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Abstract 

Physical  and  numerical  modeling  of  turbulent  mixing  and  turbulent  com¬ 
bustion  in  a  Scramjet  combustion  chamber  is  studied.  A  high-order  hnite 
volume  scheme  is  applied  to  the  flow,  accompanied  with  dynamic  subgrid- 
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scale  model  and  local  averaging  procedure  to  account  for  the  effect  of  unre¬ 
solved  small  length  scales  turbulent  fluid  behavior.  The  Arrhenius  law  kinetic 
mechanism  is  applied  to  describe  the  nonlinearity  of  chemical  reactions.  To 
model  combustion,  the  hnite  rate  chemistry  is  developed  to  characterize  the 
combustion  process  inside  the  combustion  chamber.  Results  are  compared 
among  the  hnite  rate  model,  the  hamelet  model/progress  variable  approach 
adopted  by  Stanford  PSAAP  center  and  experimental  data.  Since  the  mesh 
requirement  to  resolve  the  hame  front  is  too  strict  for  most  combustion  sim¬ 
ulation  applications,  the  thickened  hame  model  that  artihcially  expands  the 
hame  front  is  investigated.  The  ehectiveness  of  the  model  is  studied  in  a 
one-dimensional  context  and  applied  to  the  three-dimensional  Scramjet  sim¬ 
ulation.  The  feasible  level  of  thickening  factor  is  determined  to  assume  that 
the  major  hame  features  should  be  preserved.  The  model  is  further  extended 
to  dynamic  thickened  hame  to  account  for  the  coexistence  of  premixed  and 
dihusion  hame  in  the  combustion  chamber.  It  is  extended  to  a  reduced  chem¬ 
ical  mechanism  to  reduce  the  time  complexity  needed  to  solve  chemistry. 
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Chapter  1 


1  Introduction 


1.1  Background  and  Overview 

Scramjet  (Supersonic  Combustion  Ramjet)  is  a  type  of  engine  that  op¬ 
erates  under  supersonic  airflow  conditions.  The  efficiency  in  its  propulsion 
system  over  ramjet  has  made  it  a  very  active  research  topic  over  the  past 
decades.  The  flow  structure  inside  a  scramjet  combustion  chamber  is  an 
extension  of  the  jet  in  crossflow  (JICF)  configuration,  which  categorizes  a 
type  of  flow  regime  that  is  frequently  encountered  in  many  practical  appli¬ 
cations.  For  many  years,  a  large  body  of  literature  has  been  devoted  to  the 
study  of  the  structure  of  the  jet  in  crossflow  regime.  It  has  been  conhrmed 


by  experiments  36  38l  38l  41  and  numerical  studies  19,37  that  the  major 


flow  and  vortical  features  of  a  typical  jet  in  cross-flow  can  be  represented 
by  Figure  |1.1  As  shown  in  the  figure,  the  flow  structure  consists  of  an 
orihce  from  which  a  flow  is  injected  perpendicularly  into  a  crossflow.  The 
most  dominant  features  are  the  counter-rotating  vortex  pair  (CVP),  a  vortex 
structure  that  starts  forming  partly  attributable  to  the  jet  shear-layer  insta¬ 
bilities  near  the  point  of  injection  and  takes  up  much  of  the  downstream  flow 
region,  the  horseshoe  vortices  that  generates  around  the  nozzle  and  wraps 
up  the  part  of  jet  flow  just  issuing  out  from  it,  the  shear  layer  vortices  that 
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Counler-rolaling 
vortex  pair 


Figure  1.1:  The  conceptual  flow  and  vortical  features  that  can  be  identified 
in  a  jet  in  cross-flow.  Figure  adopted  from 


41 


result  from  a  Kelvin-Helmholtz  instability  near  the  nozzle  (more  extensively 


studied  in  21,36,^)  and  the  wake  vortices  that  originate  from  the  boundary 
layer  of  the  wall  |^.  Moreover,  when  the  crossflow  is  supersonic,  as  is  the 
case  in  the  Scramjet  configuration,  some  additional  structures,  including  the 
bow  shock  formed  due  to  the  blockage  of  supersonic  crossflow  by  the  injected 
inflow,  the  shock  train  structures,  and  the  recirculation  zones  generated  both 


in  the  upstream  and  downstream  vicinity  of  the  nozzle  can  be  observed  26 


The  scramjet  combustion  research  has  been  conducted  in  many  parts 
of  the  world,  both  experimentally  and  numerically.  Among  them,  ground 
and  flight  testings  that  achieve  supersonic  combustion  were  conducted  and 


post-flight  data  studied  at  the  University  of  Queensland  39  ,  a  1:1  model 
for  the  same  conditions  as  for  the  flight  experiment  was  performed  by  the 
German  Aerospace  Center  [^.  On  the  side  of  computational  methods,  35 


used  a  Reynolds-averaged  Navier-Stokes  solver  to  study  the  effect  of  different 


2 


equation  of  states  and  turbulence  models.  Large  eddy  simulation  has  also 
been  a  very  popular  approach  to  model  turbulence 


18,19.  A  number 


of  combustion  models  including  the  Eddy  Dissipation  Concept  model  11 


Assumed  Probability  Density  Function  (PDF)  approach  31,32  ,  Linear  Eddy 
Mixing  (LEM)  model  and  the  Flamele  model  has  been  proposed  to  character¬ 
ize  the  turbulent  combustion  process.  Among  them,  the  flamelet  approach 
and  its  extensions  based  on  the  assumption  that  the  flame  is  locally 
one  dimensional  are  considered  preferable  in  many  applications  due  to  its 
ability  to  decouple  chemistry  from  hydrodynamics  and  the  decrease  in  the 
computational  cost  comparing  to  traditional  models.  In  this  work,  an  LES 
solver  together  with  an  alternative  combustion  approach  called  the  hnite  rate 
chemistry  model  is  put  forward  and  discussed  in  detail.  The  advantage  of 
the  chemistry  model  here  is  that  it  does  not  rely  on  the  one  dimensional 
structure  assumption  underlying  the  flamelet  model.  With  a  simpler  struc¬ 
ture  and  implementation  effort,  it  allows  chemistry  process  to  be  evaluated 
on  a  scale  larger  than  the  smallest  scale  of  turbulence.  The  restrictions  of 
the  model  are  also  investigated. 


1.2  Dissertation  Organization 

The  Navier-Stokes  Equation  that  governs  the  fluid  dynamics  in  general 
is  introduced  in  Chapter  2,  together  with  a  brief  discussion  of  the  Reynolds 
number  as  a  quantity  that  characterizes  the  level  of  turbulence  of  a  given 
flow.  Chapter  3  goes  on  to  discuss  the  concepts  in  turbulence  and  the  various 


3 


approaches  to  numerically  modeling  turbulent  flows.  Much  focus  is  directed 
to  the  dynamic  LES  (Large  Eddy  Simulation)  approach,  which  provides  a 
good  balance  between  computational  cost  needed  and  the  level  of  resolution 
that  could  be  obtained.  Chapter  4  introduces  the  WENO  scheme,  a  high- 
order  finite  volume-based  numerical  scheme  that  is  applied  in  the  Scramjet 
numerical  simulation.  The  combustion  models  are  presented  in  Chapter  5. 
The  chemical  kinetics  central  to  the  combustion  process  is  discussed.  The 
flame  properties  and  structures  are  the  major  focus  of  a  combustion  model. 
We  discusses  two  types  of  approaches  to  characterizing  a  flame,  namely,  the 
flamelet  model  and  the  finite  rate  chemistry.  The  flame  resolution  is  studied 
in  this  chapter  in  a  one-dimensional  framework.  Chapter  6  gives  the  results 
of  the  three-dimensional  Scramjet  simulation,  followed  by  the  conclusions 
and  discussions  in  Chapter  7. 
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Chapter  2 


2  Governing  Equations 

2.1  Euler’s  Equations 

The  fundamental  governing  equations  in  fluid  dynamics  are  the  Euler’s 
equations,  which  are  derived  from  the  conservation  laws  of  mass,  momentum 
and  energy.  Euler’s  equations  characterize  flows  with  no  viscosity.  It  has  the 
differential  form 


dp  dpui 
dt  dxi 

=  0, 

(2.1) 

d  d  dp 

^tPUi  +  ^pu,u,  +  ^ 

=  0, 

(2.2) 

dpe  d  .  , 

dt 

=  0. 

(2.3) 

The  above  equations  are  written  using  Einstein  notation,  in  which  matched 
pairs  of  indices  refers  to  a  summation. 

2.2  Navier-Stokes  Equations 

It  is  assumed  in  Euler’s  equations  that  forces  exerted  on  a  surface  of  flow 
are  normal  to  that  surface.  In  other  words,  viscous  forces  are  not  taken 
into  consideration.  When  more  general  fluids  are  considered,  such  an  as¬ 
sumption  becomes  incorrect.  A  term  characterizing  viscous  forces,  the  stress 
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tensor  term  a  is  introduced,  and  the  more  general  term  —p{x,  t)n  +  a{x,t)  -n 
replaces  —p(x,  t)n  in  the  integral  form  of  conservation  of  momentum  equa¬ 
tion  [^.  With  a  few  more  assumptions  on  the  properties  of  a,  including  it 
being  linearly  dependent  on  velocity  gradient  and  being  symmetric,  a  can  be 
transformed  and  decomposed  into 

a  =  2p[D  —  -(divn)/]  -|-  ^(divu)/.  (2.4) 

o 


Here,  /  is  the  identity  matrix,  p  is  called  the  hrst  coefficient  of  viscosity,  and 
(  the  second  coefficient  of  viscosity. 

Upon  further  applying  the  divergence  theorem,  and  leaving  out  the  second 
coefficient  of  viscosity,  the  conservation  of  momentum  equation  could  take 
on  the  form 

where  Tij  characterizes  the  deviatoric  stress  tensor,  and  has  the  form 


hij  2jpSij , 


Si,  =  - 

^  2 


dui  du, 

— f  H - 2 

dxj  dxi 


1  dui 


Sii. 


3dx 


(2.6) 


The  Reynolds  number  is  a  very  important  dimensionless  quantity  in  fluid 
mechanics.  It  is  dehned  by  the  characteristic  length  of  flow  region  L,  the 

characteristic  velocity  U,  and  the  coefficient  z/  =  —  of  kinematic  viscosity, 

LU  ^ 

and  can  be  written  as  Re  =  - .  The  Reynolds  number  can  be  further 

V 

written  as  the  ratio  of  inertia  forces  to  viscous  forces,  and  is  often  used  to 
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describe  whether  a  given  flow  is  laminar,  transient  or  turbulent.  A  flow 
regime  with  a  high  Reynolds  number  {Re  >  4000)  indicates  that  viscous 
forces  are  negligible  comparing  to  inertial  forces,  and  therefore  the  flow  tends 
to  be  turbulent,  as  is  the  case  of  Scramjet  injected  flow  and  crossflow. 
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Chapter  3 


3  Turbulence  Modeling 


Turbulence  is  ubiquitous  in  nature  and  central  to  many  applications  in 
engineering,  and  is  characterized  by  chaotic,  stochastic  property  changes. 
It  is  also  well  known  as  the  most  important  unsolved  problem  of  classical 
physics. 

It  is  discovered  that  the  details  of  turbulence  itself  appear  to  be  un¬ 
predictable,  but  its  statistical  properties  are  reproducible.  Although  the 
Navier-Stokes  equations  system  is  deterministic,  a  statistical  description  of 
turbulence  is  needed,  especially  for  high-Reynolds  number  flow,  due  to  the 
fact  that  fluid  flows  are  very  sensitive  to  changes  in  initial  and  boundary 
conditions,  and  can  become  unstable  and  irregular  in  time  and  in  space. 

The  Russian  mathematician  A.  N.  Kolmogorov  contributed  a  great  deal  in 


the  development  of  turbulence  research.  In  his  1941  work  22  ,  later  referred 


to  as  K41  theory,  he  made  several  hypotheses  which  laid  down  the  frame¬ 
work  of  turbulence  study  afterwards.  One  of  them  is  that  when  Reynolds 
number  for  a  fluid  flow  is  very  high,  small  scale  turbulent  motions  become 
isotropic.  This  is  not  true  for  flows  of  large  scales,  because  the  flow  motions 
are  associated  with  geometrical  features.  Isotropic  small  scales  are  the  same 
for  all  turbulent  flows.  In  addition,  he  postulated  that  this  small  scale  is  a 
function  only  of  molecular  viscosity  v  and  energy  dissipation  rate  (per  unit 


mass)  e.  Since  the  quantity  v  has  dimension  [-h]^[T]  where  [L\  and  [T]  are 
units  of  length  and  time,  respectively,  and  e  has  dimension  [L]^[T]“^,  using 
dimensional  analysis,  the  scale  rj  can  be  represented  hy  7]  = 

His  other  hypothesis  is  that  at  very  high  Reynolds  numbers,  all  the  scale 
statistical  properties  between  rj  and  L,  the  characteristic  length  of  flow,  can 
be  universally  and  uniquely  determined  by  r  (the  scale),  and  the  energy 
dissipation  rate  e. 

Here,  r]  is  the  so-called  Kolmogorov  length  scale.  For  large  scale  vor¬ 
tices,  viscous  effects  are  not  important,  and  kinetic  energy  of  turbulence 
is  mostly  contained  in  those  structures,  and  transferred  from  larger  scales 
to  smaller  scales.  When  scales  are  small  enough  to  reach  the  Kolmogorov 
length  scale,  molecular  diffusion  becomes  important  and  viscous  dissipation 
of  energy  takes  place,  turning  kinetic  energy  into  internal  energy.  In  other 
words,  in  the  range  between  characteristic  length  scale  L  of  the  flow  and 
Kolmogorov  length  scale  rj,  inertial  effects  play  a  much  more  important  role 
than  viscous  effects,  and  it  makes  sense  to  assume  that  viscosity  does  not 
play  a  role  within  it.  This  range  is  called  inertial  range. 

The  way  kinetic  energy  is  distributed  among  different  scale  of  turbulent 
flow  can  be  characterized  by  the  energy  spectrum  function  E{k)  which  sat- 
ishes 

E{k)dk  =  total  kinetic  energy,  (3.1) 

where  k  is  the  modulus  of  the  wavevector  k  corresponding  to  the  transfor- 
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Figure  3.1:  Energy  spectrum  in  the  time  domain  for  data  from  SI.  Reynolds 
number  R\  =  2720.  Courtesy  Y.  Gagne  and  M.  Marchand 


mation  of  quantities  from  the  physical  space  to  the  Fourier  space,  and  the 
wavenumber  fc  of  a  vortex  of  spatial  dimension  r  is  given  hy  k  =  There¬ 
fore,  E{k)dk  represents  the  contribution  to  total  kinetic  energy  of  the  Fourier 
harmonics  having  the  absolute  value  of  their  frequency  between  k  and  k  +  dk. 
Using  Kolmogorov’s  third  hypothesis  and  dimensional  analysis,  in  the  inertial 
range, 

E{k)  =  (3.2) 


This  is  the  famous  —  |  law 
most  famous  experimental 


for  the  energy  spectrum.  Figure 


3.1 


is  one  of  the 


evidence  that  support  this  law,  taken  from  14 
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3.1  Reynolds- Averaged  Navier-Stokes  Equations 

RANS  equations  are  time-averaged  equations  typically  used  to  model 
turbulent  flows.  They  involve  decomposition  of  instantaneous  variables  into 
its  mean  and  fluctuation  part:  u  =  u  +  u' .  After  imposing  the  time-averaged 
Reynolds  operator  to  the  governing  equations,  the  Navier-Stokes  equations 
can  take  the  form 


d{fmi) 

dt 

djpE) 

dt 


+ 

+ 


dp 

m 

d 

dxj 


+ 


dpuj 

dxj 


(yPUiUj  ) 


A 

dxj 


{ujpH) 


0, 

dp  d^  _  d{pTij) 

dxi  dxj  dxj 


(3.3) 

(3.4) 

(3.5) 


Here,  u  is  the  Favre-averaging  defined  for  any  variable  u  such  that  u  =  —, 

P 

and  a  corresponding  decomposition  gives  u  =  u  +  u'  =  u  +  u” .  Besides,  in 


11 


the  above  equations,  Tij  =  u"u",  and 


a, 


E 

H 


(Ij 


pE”u" 


/  1  \ 

(3.6) 

2/i  (  Sij  1 , 

^  ,  UiUi  , 

"  2  2  ’ 

(3.7) 

^  +  -, 
p 

(3.8) 

dT 

(3.9) 

dxj 

Cppu''T"  +  Ui{pTij  -  aij) 

pu'Ju'Ju" 

+  2 

Terms  that  includes  the  nonlinear  Reynolds  stress  term  pvij  requires  fur¬ 
ther  modeling  efforts  in  order  to  bring  a  closure  to  the  RANS  equations. 


3.2  Direct  Numerical  Simulation 

Direct  numerical  simulation  (DNS)  is  an  approach  used  to  solve  Navier- 
Stokes  equation  without  resorting  to  any  turbulence  models.  Application  of 
DNS  requires  that  all  temporal  and  spatial  scales,  from  the  characteristic 
length  scale  of  the  flow,  down  to  the  Kolmogorov  scales,  are  to  be  resolved. 
The  choice  of  numerical  schemes  is  essential  to  the  accuracy  of  DNS,  thus  low 
dissipation  and  dispersion  error  schemes  are  to  be  sought.  Though  DNS  is 
considered  to  give  very  detailed  and  accurate  results,  the  computational  cost 
to  meet  the  corresponding  resolution  requirements  can  be  prohibitive.  The 
cost  becomes  even  more  expensive  when  flows  with  high  Reynolds  number 
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is  to  be  simulated.  It  is  estimated  that  for  three-dimensional  simulations, 
the  number  of  operations  could  be  proportional  to  as  high  as  Re^.  Because 
of  this,  DNS  is  considered  more  as  a  research  tool  running  on  small-scale 
problems,  and  not  suitable  for  solving  practical  engineering  problems. 

3.3  Large-Eddy  Simulation 

Large-Eddy  Simulation  (LES)  is  another  class  of  popular  approaches  for 
simulating  turbulent  flows.  Unlike  RANS,  the  flow  held  is  resolved  down  to  a 
certain  length  scale,  and  solved  via  hltered  Navier-Stokes  Equations.  Scales 
smaller  than  that  are  called  subgrid  scales.  These  scales  are  eliminated  from 
the  governing  equations,  and  they  are  modeled  rather  than  resolved.  The 
resolved  scales  are  usually  much  larger  than  the  Kolmogorov  scale,  therefore 
making  LES  much  more  manageable  computationally  than  DNS,  and  gener¬ 
ally  has  a  complexity  in  between  DNS  and  RANS.  LES  involves  separating 
velocity  held  into  resolved  and  sub-grid  components.  The  resolved  part  is  a 
hltered  held,  and  the  hltering  operation  can  be  dehned  generally  in  terms  of 
a  convolution  against  a  kernel  F 


u  =  F  *  u  (3.11) 

where  n  is  a  how  held  variable.  Among  diherent  hltering  operations,  a 
density-weighted  Favre  hltering  is  widely  used,  with  the  purpose  of  avoid¬ 
ing  modeling  the  mass  concentration  equation  and  simplify  the  resulting 
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equation  system.  For  an  arbitrary  quantity  q,  Favre  filtering  is  defined  as 

~  pq 
q  =  —■ 

P 

After  applying  this  filter,  the  Navier-Stokes  equations  have  the  following 
form 


^  dpuj 
dt  dxi 
dpui  ^  dpuiUj 


dt 


dxi 


dt 


dxj 


=  0, 


dE  d{E  +  p)uj 
dt  dxj 


dpnp  ^  dpnpUj 


ddijUi  d 
dxj  ^ 


Hi)pD^ 


dxj 


dxj  ’ 


where  the  Subgrid  Scale  terms  are  expressed  as 


Tij 

(3.16) 

(H) 

q) 

=  p{cpTuj  -  Cpfuj), 

(3.17) 

(T) 

q) 

(3.18) 

o'r 

(3.19) 

if' 

fi 

1 

(3.20) 
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and  the  filtered  variables  are  expressed  as 


E 

,  pVk  ,  Tkk 

=  pe+  2  +  2  ’ 

(3.21) 

Hk 

~  ,  P 

—  Gh  + 

P 

(3.22) 

Hi 

~  ,  P 

=  ei  + 

(3.23) 

P 


Here,  E,  Hh  and  Hi  stand  for  total  specific  energy  and  partial  specific  en¬ 
thalpy  of  two  species,  respectively. 

3.3.1  Dynamic  Smagorinsky  Model 

After  the  equations  are  modihed  and  variables  all  transformed  into  hltered 
quantities,  the  next  step  is  to  model  those  SGS  terms  left  in  these  equations. 
The  term  Tij  that  appears  in  momentum  equations,  which  is  also  known  as 
SGS  stress  tensor  is  the  hrst  crucial  term,  and  the  only  one  that  needs  to  be 
modeled  in  order  to  close  the  momentum  equation. 

This  term  can  be  separated  into  a  trace-free  term  and  an  isotropic  term 

Tij  =  r“.  -f  tP  =  (^Tij  -  TkkY^  +  Tkk^-  (3-24) 

For  incompressible  flows,  only  the  anisotropic  tensor  is  relevant  since  the 
isotropic  part  is  absorbed  into  the  pressure;  for  compressible  flows  on  the 
other  hand,  both  terms  should  be  considered.  With  the  trace-free  Smagorin- 
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sky  eddy  viscosity  model, 


r-,  =  n,  -  rj-f  =  -2CspA^\S\{S,,  -  (3.25) 

where  Sij  =  ^  1*^1  ^  ^ 

filter  width.  With  Yoshizawa’s  expression, 

rlj  =  TkkY  =  ‘^SijCiA‘^p\S\^.  (3.26) 

Therefore,  the  modeled  subgrid  scale  (SGS)  stress  Tij  is 

=  -2Cs-pA^\S\(S,^  -  \SkA,)  +  ^-5,,CiA^-p\S\\  (3.27) 

In  order  to  close  the  momentum  equation,  hrstly,  using  Favre-hltering 
relation,  Tij  can  be  further  written  as 

_  _ _  _  pUi  pUj 

Tij  =  p[UiUj  —  UiUj)  =  pUiUj - ^ (3.28) 

In  the  dynamic  SGS  model,  the  central  idea  is  the  introducion  of  a  test  hlter 
which  is  larger  than  the  resolved  grid  hlter  (in  current  Scramjet  problem,  it 
is  set  to  be  two  times  the  dimension  of  the  original  grid  hlter,  but  in  fact  can 
be  chosen  at  will)  so  that  the  data  information  can  be  utilized.  The  hlter 
width  of  the  test  hlter  is  denoted  by  A.  The  test  hlter  is  always  applied  to 
the  resolved  held,  and  on  this  hlter,  a  stress  term  similar  to  the  form  of  r^ 
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can  be  defined  as  follows 


rp  _ -  pUipUj 

p 


(3.29) 


With  Germano’s  Identity  24  ,  the  Leonard  stress  tensor  Lj,-  is  expressed 


in  terms  of  and  Tij 


Lij  —  Tij 


T, 


pUiVjj 


pUkPUj 

"p 


(3.30) 


/  pUj  pUj  \  _  pUkpUj 

VP/  P 


pUipVjj 


(3.31) 

(3.32) 


The  Leonard  stress  is  totally  resolved,  while  the  other  two  stress  tensors  in 
this  relation,  and  r^-  are  both  modeled  using  the  trace-free  Smagorinsky 
eddy  viscosity  model,  and  we  obtain  the  following  relation 


L%  =  T-  -  =  2CsA^p\S\S^^  -  2CsA^^\S\Si^  =  CsM^p 


(3.33) 


where 

=  2A‘^ p\S\St^  -  2A^m%.  (3.34) 

The  above  equation  is  actually  an  overdetermined  system  of  equations 
which  corresponds  to  five  (off-diagonal  components  of  stress  matrix)  relations 
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between  M“-  and  L“-.  In  order  to  give  an  approximate  solution  for  Cs,  the 
most  common  way  is  via  least-square  method.  When  a  target  region  is  chosen 
(details  explained  below),  an  average  over  variables  and  within  this 
region  is  applied  and  the  resulting  formula  for  Cs  becomes 


((E 14, M^)*) 

((EA/gAf§)+)' 


(3.35) 


Here,  the  expression  means  only  positive  terms  are  added  up,  and  neg¬ 
ative  terms  are  ignored.  The  expression  (...)  indicates  the  averaging  process. 

In  order  to  determine  coefficient  Cj,  a  similar  process  leads  to  the  expres¬ 
sion 


T  I  _  7^2 

'^kk  kk 


fl,  =  -2CjA^p\S\^  +  2CiA^p\S\^  =  CjMl,  (3.36) 


Ml,  =  -2A^p\S\^  +  2A^p\S\^. 


(3.37) 


Similar  to  the  above  treatment,  since  this  relation  corresponds  to  three 
(diagonal  components  of  stress  matrix)  equations,  an  averaging  over  a  spec¬ 
ified  area  leads  to  the  least-square  type  of  solution  of  Cj  in  the  form  of 


Ci  = 


(ML) ' 


(3.38) 


( M)  ('V') 

In  the  energy  equation,  among  the  three  sub-grid  terms  g)  ,  g)  and 
the  heat  transport  flux  gj'^^  is  the  one  with  the  lowest  order,  and  the 
only  one  considered  here  in  the  modeling  process.  (  gj^^  and  gj^^  are  set  to 
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be  zero).  Heat  transport  flux  is  modeled  as 


(H)M 

Qi  =  -pCp 


df  CpUt  df 


Prt  dxi 


Prt  dxi  dxi 


dT 

— ,  (3.39) 


where  Prt  is  the  Prandtl  number  and  the  one  to  be  determined  by  the 
dynamic  model  in  a  similar  fashion. 

Applying  Favre  filtering  to  the  flux  term  gives 


=  pc  Tv  - 

p 


(3.40) 


and  is  modeled  by 


riH)M  ,.CsA2|S|ar 


(3.41) 


Germano’s  identity  again  gives  rise  to  the  following  relation 


Lf^  =  gf _pcp^m  ^ (342) 


p- 


p- 


Eddy  diffusivity  SGS  model  for  Pt  gives 


(3.43) 


where 


Mf'  =  AMp4IS|f  -AyaAlf. 


(3.44) 


Again,  a  least-square  approach  is  applied  to  the  system  of  three  equations 
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and  the  Prt  is  calculated  as 


Pr.  =  C< 


(3.45) 


Lastly,  we  model  the  term  in  species  concentration  equations  and 
derive  the  Schmidt  number  Sct  dynamically.  The  linear  gradient  assumption 
leads  to 

(3-46) 

At  test  hlter  grid  level, 


=  pvi-ip  - 


pvipjj 

P 


(3.47) 


is  modeled  by 


Qr  =  -p 


Sct  dxj 


(3.48) 


Using  Germano’s  identity  again 


P 


P 


P 


(3.49) 


The  modeled  quantity  is 


where 


(3.51) 


The  modeled  Schmidt  number  Set  is  thus  calculated  via  least-square  ap¬ 
proach 

Set  =  h,  .  (3.52) 

((Eh*’v"’)+) 


3.3.2  Averaging  Procedure 

Before  least-square  approach,  the  dynamic  SGS  model  used  in  the  pre¬ 
vious  section  to  derive  the  expression  of  the  eddy  viscosity  coefficients,  the 
Prandtl  number  and  the  Schmidt  number  are  all  functions  of  time  and  space. 
In  various  applications,  however,  it  is  shown  that  the  the  right-hand  side  of 
the  htting  equation  could  become  zero  in  some  cells.  In  order  for  the  model 
to  be  consistent  and  well  conditioned,  the  coefficients  in  question  are  com¬ 
monly  averaged  over  a  specihed  region  in  the  flow.  The  way  to  determine 
the  appropriate  area  of  averaging  is  problem  specihe,  and  is  usually  depen¬ 
dent  on  the  properties  of  the  flow.  The  concept  of  statistical  homogeneity 
of  flow  is  introduced,  and  it  is  assumed  that  coefficients  Cs,  Cj,  Pr  and  Sc 


are  independent  in  the  homogeneous  directions  of  the  flow  28  .  Therefore, 


an  averaging  procedure  can  be  applied  along  the  homogenous  direction  to 
adjust  the  values  of  coefficients  away  from  zero.  In  the  context  of  Rayleigh- 
Taylor  instability  problems,  for  example,  a  certain  level  of  homogeneity  exists 
on  all  points  in  space  with  the  same  z-coordinate,  and  therefore,  planar  av- 
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Figure  3.2:  H2  concentration  profile  outside  of  nozzle  for  240  micron  (left) 
and  120  micron  (right),  using  local  averaging  procedure. 


erage  could  be  applied  on  the  dynamically  calculated  coefficients.  For  the 
3D  Richtmyer-Meshkov  instability  problem  with  a  circular  geometry,  on  the 
other  hand,  a  spherical  average  is  usually  employed.  In  the  case  of  Scramjet 
simulation,  however,  the  complexities  of  flow  on  many  scales  means  that  a 
homogeneous  direction  may  not  exist  after  all.  Therefore,  the  model  is  fur¬ 
ther  adapted  so  that  the  region  upon  which  to  apply  averaging  is  restricted 
to  small  mesh  blocks  whose  dimensions  are  left  to  choice,  called  localized 
mesh  block  averaging.  The  advantage  of  this  adaptive  SGS  model  is  that 
it  no  longer  involves  data  transfer  and  communications  among  processors, 
which  could  be  a  bottleneck  for  either  planar  or  spherical  averaging.  The 
freedom  to  choose  the  size  of  averaging  mesh  blocks  also  makes  the  model 


more  flexible  and  attractive  to  use.  Figure  3^  illustrates  the  effect  of  local 
averaging  procedure  on  the  Scramjet  simulation  with  two  grid  levels. 


22 


3.3.3  Traditional  Smagorinsky  Model 


Comparing  to  the  dynamic  version  of  the  sub-grid  model  described  in  the 
above  section,  the  traditional  Smagorinsky  model  assumes  the  eddy  viscosity 
coefficients  are  constant,  and  should  be  chosen  differently  for  different  flow 
helds.  Moreover,  there  are  some  other  restrictions.  According  to  [^,  (i) 
the  model  does  not  have  the  correct  limiting  behavior  near  the  wall;  (ii) 
the  model  does  not  vanish  in  laminar  flow,  and  it  is  demonstrated  to  be 
too  dissipative  in  the  laminar-turbulent  transition  region;  (hi)  the  model 
does  not  account  for  backscatter  of  energy  from  small  scales  to  large  scales, 
which  has  been  shown  to  be  of  importance  in  the  transition  regime;  (iv) 
compressibility  effects  are  not  included  in  the  model.  Dynamic  SGS  model, 
based  on  the  application  of  Smagorinsky  model  at  two  different  hlter  level,  is 
therefore  developed  to  overcome  some  of  those  drawback,  and  produce  more 
accurate  simulation  results.  However,  the  Smagorinsky  model  is  still 
widely  applied  in  many  large  eddy  simulation  problems  due  to  its  simplicity 
and  effectiveness  in  many  applications. 
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Chapter  4 


4  Numerical  Methods 


4.1  Finite  Volume  Approach 

The  finite  volume  approach  is  a  numerical  method  based  on  the  integral 
form  of  a  PDE.  It  is  well  suited  as  a  numerical  scheme  to  describe  physical 
properties  of  flows.  Conservation  laws  are  a  very  important  class  of  homoge¬ 
neous  hyperbolic  partial  differential  equations,  and  they  are  usually  derived 
directly  from  physical  laws.  A  one-dimensional  conservation  law  has  the  form 


Ut(x,t)  +  f(U(x,t)),  =  0. 


(4.1) 


Here,  f{U{x,t))  is  called  the  flux  function,  and  U  is  an  unknown  function. 


usually  a  certain  flow  quantity.  Let  Xj  =  jAx,  Equation  4.1  can  be  further 


written  into  an  integral  form  over  an  interval  [xj  —  ^Ax,Xj  +  \Ax\ 


^-L  ‘A  +  Iai,; 


(4.2) 


1  CC  *  I  — 

Let  Qj  =  ^  f{U{x,t))dx,  and  averaging  over  +  At]  gives 


Q 


U+At  _  U 

At 


1  rti+At  ^  ^  nti+At 

=  Xti  nU{xi-^Ax.t))dt-^J^^ 


f{U{xj+-Ax,t))dt. 

(4.3) 
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(4.4) 


=  f{U{xj  —  ^Ax,t))dt,  after  rearranging 


)ti 


ti+At 


If  we  can  approximate  the  average  flux  F*_i  and  ,,  we  obtain  a  discrete 
numerical  scheme.  The  choice  of  this  numerical  flux  function  dehnes  various 
forms  of  hnite  volume  schemes. 


4.2  WENO  Scheme 

The  WENO  scheme  is  a  high-resolution  hnite  volume  scheme  developed 
to  combine  stencils  of  lower  orders  to  reach  a  higher  order  approximation. 
It  extends  the  ENO  scheme  which  selects  a  single  candidate  stencil  only. 
WENO  scheme  uses  a  convex  combination  of  all  candidate  stencils,  and  each 
one  is  given  a  weight  that  corresponds  to  the  contribution  of  that  stencil  to 
the  total  numerical  hux.  The  weights  are  dehned  so  that  away  from  how 
discontinuities,  it  approaches  certain  optimal  weights  in  order  to  get  a  higher 
order  accuracy;  in  regions  of  discontinuities  and  shocks  on  the  other  hand, 
the  stencils  that  include  the  discontinuous  points  are  assigned  a  near  zero 


weight  20 


In  the  following,  a  5*^  order  hux  approximation  of  the  WENO  scheme 
is  described,  largely  following  the  derivation  of  [^.  The  scheme  is  the  one 
being  used  in  Scramjet  simulation  to  achieve  a  high  order  approximation  of 
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complicated  flow  structures.  The  flux  approximation  can  be  expressed  as 

2 

f  =  '^Wkhk{f-2,f-i,fo),  (4.5) 

where  hk  are  simply  linear  combinations  of  fluxes  within  the  stencil  k,  and 
can  be  expressed  as  hk{f-2,  f-i,  fo)  = 

For  stability  purposes,  the  approach  described  below  are  applied  to  /■*■ 
and  f~  respectively,  where  f~^  and  f~  are  the  positive  and  negative  compo¬ 
nents  corresponding  to  flux  splitting  f{x)  =  f~^{x)  +  f~{x),  where  —f~^{x)  > 

(JjJb 

0,  and  <  0.  The  Lax-Friedrichs  flux  splitting  is  the  one  commonly 

ax 

used,  and  has  the  following  form 

f^{x)  =  ^{f{x)±ax),  (4.6) 

where  a  is  chosen  to  be  the  global  maximnm  of  |/'(a;)|. 

For  a  system  of  eqnations,  the  WENO  scheme  and  the  corresponding 
weights  of  candidate  stencils  are  compnted  in  the  local  characteristic  fields. 
The  Jacobian  matrices  are  evaluated  at  the  average  mean  values  of  state  vari¬ 
ables  of  two  neighboring  cells.  Denote  h  and  r*  as  the  left  and  right  eigen¬ 
vectors  of  the  Jacobian  matrix  Jj  at  local  characteristic  held  i,  respectively. 
The  approximation  of  flux  in  the  mapped  held  can  therefore  be  replaced  by 

2 

fi  =  '^Wkhk,  (4.7) 

A:=0 
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r 

k 

1  =  0 

1  =  1 

1  =  2 

3 

0 

1/3 

-7/6 

11/6 

1 

-1/6 

5/6 

1/3 

2 

1/3 

5/6 

-1/6 

Cl 

k  =  0 

k  =  1 

k  =  2 

r  =  2 

r  =  3 

1/3 

1/10 

2/3 

6/10 

3/10 

Table  4.1:  Coefficients  a^  i  (left)  and  Cl  (right) 

where  Wk  =  Wfc(h-f-2,  h'f-i,  h-fo,  h'fi,  li-f2)  and  hk  =  hfc(h-ffc-2,  h'ffc-i,  h'ffc) 
at  stencil  k.  Wk  is  dehned  by 


Wk 


o^k 

Qfo  +  Q?!  +  Cti 


(4.8) 


where 


Cl 


Oik  — 


{e  +  ISk)P^ 


k  =  0,1,2 


and  ISk-,  {k  =  0, 1,  2)  are  evaluated  to  be 


(4.9) 


ISo  = 

1 0  1 

^(/-2  —  2/_i  +  /o)^  +  ^(/-2  —  4/_i  +  3/0)^, 

(4.10) 

ISi  = 

1 3  1 

-(/_!-  2/0  +  /l)^  +  -(/_l-/l)^ 

(4.11) 

IS2  = 

1 3  1 

-(/o-2/i  +  /2)2  +  -(3/o-4/i  +  /2)1 

(4.12) 

Cl  is  given  in  Table  4.1,  e  is  a  small  positive  value  introduced  to  avoid  the 


denominators  of  being  zero.  The  power  p  is  chosen  to  be  2  as  20  suggests. 
Fluxes  /  are  then  mapped  back  to  the  spatial  held  by 


/  =  (4.13) 

i=l 
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using  right  eigenvectors  of  Ji. 


The  time  discretization  L{U)  = 


qr\U)-q^^iU)  , 


At 


is  approximated  by  a 


3th  TVD  Runge-Kutta  scheme  expressed  by 


=  U^  +  AtL{U^), 

(4.14) 

=  lu" +  \u'" +  \l^tL(U\ 

(4.15) 

Un+l 

=  T"  + + -AiL(t;P>). 

3  3  3 

(4.16) 

Chapter  5 


5  Combustion  Models  and  Mesh  Resolution 


Analysis 


5.1  Thermodynamic  Properties  and  Equations  of  State 

The  interactions  among  reacting  species  in  both  hydrodynamic  and  com¬ 
bustion  process  call  for  an  accurate  description  of  their  thermodynamic  prop¬ 
erties  and  multi  species  equations  of  state.  Most  formulas  and  results  pre¬ 
sented  in  this  section  are  widely  used  in  the  thermodynamics  community, 
and  can  be  found  in  standard  references  ,  For  simulation  purposes, 

these  equations  of  state  and  thermodynamic  relations  are  implemented  and 
grouped  into  a  library  separate  from  the  hydrodynamic  part  of  scheme  to 
allow  for  fast  access  from  the  hydro  code. 


The  theory  of  thermodynamics  33  indicates  that  for  any  given  medium. 


only  two  of  the  parameters  P,  T,  V,  E  and  S  are  independent.  A  fundamental 
relation  in  thermodynamic  theory  (derived  from  the  hrst  and  second  law  in 
thermodynamics)  involves  the  change  in  internal  energy.  To  be  more  specihc, 
a  change  in  the  internal  energy  in  a  portion  of  medium  is  the  sum  of  the  work 
done  on  the  medium  in  the  form  of  pressure,  plus  the  heat  transferred  to  the 
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medium.  This  fact  is  often  expressed  as 


dE  =  TdS  -  PdV,  P{S,  V) 


dE 

W 


T{S,  V) 


dE 


V 


Another  fundamental  thermodynamic  property  of  a  flow  is  that,  given  con¬ 
stant  entropy,  the  pressure  increases  with  increasing  density,  or  in  functional 

form,  P  =  P{p,  S).  This  relation  leads  to  a  convenient  definition  of  a  positive 
dP 

quantity  (?  =  — —  =  iA(p,  S).  The  quantity  c  is  the  well  known  sound  speed. 
dp 

The  equation  of  state  is  a  set  of  equations  that  describes  the  relationship 
between  thermodynamic  state  variables  nnder  given  physical  conditions.  Two 
very  common  and  closely  related  classes  of  equation  of  state  for  a  single 
species  are  the  ideal  gas  and  a  one  parameter  family  of  ideal  gases,  the 
gamma-law  gas. 

An  ideal  gas  is  a  medinm  that  satishes  the  law  of  Boyle  and  Gay-Lussac, 
and  can  be  written  in  the  form 


PV 


RT 

’ 


(5.1) 


where  R  is  the  universal  gas  constant,  and  M  is  the  molecular  weight  of  the 
gas.  A  very  important  feature  of  an  ideal  gas  is  that  the  internal  energy 
is  a  function  of  temperature  only,  or  E  =  E{T).  A  gamma-law  gas  further 
assumes  that  the  internal  energy  is  simply  proportional  to  the  temperature, 
and  can  be  written  as  E  =  CyT,  where  constant  Cy  is  the  specihc  heat  at 
constant  volume.  When  a  gas  is  assumed  to  be  a  gamma-law  gas,  further 
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derivation  shows  that 


P  =  P(p,S)=A(S)p-'. 


(5.2) 


Here,  A{S)  is  a  function  of  entropy  alone,  and  7  is  called  adiabatic  exponent, 

and  is  a  constant  commonly  lying  within  the  range  1  <  7  <  |.  Moderate- 

temperature  air,  which  is  a  mixture  of  gases,  has  a  value  of  about  7  =  1.4. 

Further  derivation  shows  that  7  can  be  written  as  a  ratio  of  two  quantities, 

the  specihc  heat  at  constant  pressure  Cp  and  the  specihc  heat  at  constant 
Cp 

volume  Cv-  7  =  — •  The  two  quantities  are  dehned  by 


Cy 


(5.3) 


and  related  by  the  identity  Cp  —  Cy  =  R/M. 

In  the  Scramjet  simulations,  the  assumption  of  an  ideal  gas  is  adopted, 
and  the  quantities  Cp,  H,  S  can  be  written  as  functions  of  T,  so  that  Cp  = 
Cp(T),  H  =  H{T)  and  S  =  S{T).  For  simulation  purposes,  polynomial  ap¬ 
proximations  are  derived  for  thermodynamic  properties  and  transport  coeffi¬ 
cients  of  a  great  number  of  gaseous  species,  and  least-square  htting  to  exper¬ 
imental  measurements  are  used  to  calculate  polynomial  coefficients.  NASA 

provides  tables  of  these  coefficients  that  hts  a  fourth-order  polynomial  to 

Cp(T) 

the  dimensionless  quantity  — for  individual  species,  and  further  derives 


R 
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functional  forms  of  H{T)  and  S{T)  via  integration  as  follows 


Cp 

R 

—  CLi  (I2T  T  P  ci/{r^  P 

(5.4) 

H 

T  T2  06 

W 

—  aiP  02—  -|- 

(5.5) 

S 

rp2  rp'i  rpi 

R 

=  OilnT  -|-  (I2T  P  — 1“  *^^~3 — ^  ^ 

(5.6) 

To  improve  the  table  accuracy,  the  temperature  range  is  divided  into  the 
lower  part,  from  300K  to  lOOOK,  and  the  higher  part,  from  lOOOK  to  5000K. 
The  htting  requires  that  the  two  sets  of  coefficients  have  the  common  values 
at  the  joining  point  T  =  lOOOK. 

In  Scramjet  numerical  modeling,  the  internal  energy  E  and  the  density 
p  are  known  quantities,  and  the  pressure  P  is  needed  to  close  the  hydrody¬ 
namic  equations.  The  temperature  T  should  be  determined  for  the  chemistry 
calculation.  This  is  done  using  a  nonlinear  equation  solver,  which  is  derived 
using  the  expressions 

P  PT 

P  =  (7(T)  -  l)pE,  -  =  (5-7) 

7  Cp-  Cv  =  T.  (5.8) 

Therefore, 
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When  written  explicitly  in  terms  of  T 


RT 

~w 


di  +  02^  +  cb^T'^  +  +  cisT^ 

oi  +  a2T  +  osT^  +  a4T3  +  —  1/M 


-1]E 


(5.10) 


Further  simplihcations  give 


(ai  +  a2T  +  —  1/M)RT  =  E. 


(5.11) 


The  last  equation  is  solved  with  respect  to  T  with  root-hnding  algorithms. 
For  a  single  species,  using  kinetic  theory,  the  individual  /Xj  and  Aj  are 


modeled  as  16 , 25 


/ifc  =  2.6693  X  10 


5 ^/W 


.  .  „  5 

At-Ai|Cp»+-^ 


(6.12) 

(5.13) 


(2  2) 

Q/'  is  called  the  collision  integral.  It  is  a  function  of  the  reduced  tempera¬ 
ture  Tr,  defined  as  the  ratio  of  the  actual  temperature  divided  by  the  critical 
temperature,  and  is  confirmed  by  experimental  data  to  have  the  following 
empirical  expression 

^  exp[0.45667-  0.53955  ln(T^)  + 0.18265(ln(T^))2- 

0.03629(ln(T^))3  +  0.00241(ln(T^))^].  (1  <  <  90) (5.14) 


For  a  mixture  of  N  ideal  gases,  the  total  pressure  is  the  sum  of  partial 
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pressures  of  individual  gases  according  to  Dalton’s  law 


P  =  in  which  pfc  =  pfc—T.  (5.15) 

Here,  pk,  pk  and  Mk  stand  for  partial  pressure,  partial  density,  and  molecular 
weight  of  individual  gases.  Since  the  total  density  is  the  sum  of  the  partial 
density  of  the  individual  species, 

N 

P  =  '^Pk  (5.16) 

k=l 

the  average  molecular  weight  of  the  mixture  can  be  evaluated  as 


M  = 


Z^fc=i  A4 


(6.17) 


where  Yk  =  — -  stands  for  mass  fraction  of  species  /c,  ruk  is  the  mass  of  species 
m 

k  in  a  given  volume,  and  m  is  the  total  mass  in  that  volume.  A  mixture  of  7 


law  gases  is  still  a  7  law  gas.  According  to  the  derivation  in  13  ,  the  mixture 


7(T)  can  be  expressed  as  a  combination  of  individual  gamma-law  gases 


N 


mY, 


Yt 


^(7.(r)-l)A4' 


(5.18) 


where  7j(T)  is  the  value  of  gamma  of  species  i. 

The  dynamic  viscosity  and  thermal  conductivity  of  the  gas  mixture  are 
calculated  empirically  based  on  averages  of  the  individual  species  properties 
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0 


(5.19) 


Here 


and 


k=i  k=i  ^ 


N 


M. 


^k  =  J2Gkij^Yi, 


k=l 


(5.20) 


1 


( 


Mk 


-1/2 


fii)  \Mk 


(5.21) 


5.2  Chemical  Kinetic  Mechanism 

Detailed  chemistry  is  the  starting  point  for  the  calcnlation  of  combnstion 
processes.  The  detailed  chemistry  for  H2  combustion  includes  21  reversible 


elementary  reactions  30  .  This  mechanism  is  widely  accepted  and  found 


to  be  in  good  agreement  with  experiments  over  a  wide  range  of  pressure, 


temperature  and  equivalence  ratios  in  most  applications.  Table  |5.1|  shows 
the  full  21  reaction  mechanism  and  the  corresponding  rate  parameters. 

An  N  species  chemical  reaction  has  the  form 


N 


N 


k=l  k=l 


(5.22) 


where  Sk  stands  for  the  /cth  species  in  a  mixture,  and  z/(,,  v'l  denote  the  molar 
stoichiometric  coefficients  of  species  k  as  reactant  and  product,  respectively. 


35 


Table  5.1;  Detailed  chemical  kinetic  mechanism  table.  Reactions  that  include 
the  symbol  M  involve  third-body  reactants.  The  third  body  is  an  inert 
molecule  that  is  essential  for  a  specihc  reaction  to  proceed.  The  reaction  rates 
depend  on  the  concentrations  of  third  body  with  their  specihc  efficiencies. 
A,  n  and  E  are  related  hj  k  =  ,  units  are  mol/cm^,  s“^,  K  and 

KJ/mol,  respectively.  “M  =  2.50[H2\  +  12AQ[H20]  +  lA[Other].  = 

2.00[H2\  +  Q.O[H20]  +  1.0[Other] 


Reaction 

A 

n 

E 

Hydrogen-oxygen  chain 

I.H  +  O2  ^OH  +  0 

3.52  X  10^® 

-0.7 

71.4 

2.H2  +  O^OH  +  H 

5.06  X  10^ 

2.7 

26.3 

3.H2  +  OH  ^  H2O  +  H 

1.17  X  10^ 

1.3 

15.2 

A.H2O  +  0  ^  OH  +  OH 

7.60  X  10° 

3.8 

53.4 

Direct  recombination 

+  H  +  M^H2  +  M 

1.30  X  lO^^ 

-1.0 

0.0 

+  OH  +  M  ^  H2O  +  M 

4  X  10^^ 

-2.0 

0.0 

T.O  +  0  +  M^02  +  M 

6.17  X  10^® 

-0.5 

0.0 

+  0  +  M  ^  Off  +  M 

4.71  X  10^® 

-1.0 

0.0 

9“.0  +  Ohf  +  M  ^  HO2  +  M 

8.00  X  10^® 

0.0 

0.0 

Hydroperoxyl  reactions 

lO^.TT  +  O2  +  M  ^  HO2  +  M 

ko  5.75  X  lO^*^ 

-1.4 

0.0 

koo  4.65  X  10^2 

0.4 

0.0 

II.HO2  +  H  ^OH  +  OH 

7.08  X  10^® 

0.0 

1.2 

I2.HO2  -\-  H  — )■  H2  -\-  O2 

1.66  X  10^® 

0.0 

3.4 

I3.HO2  +  H  H2O  +  0 

3.10  X  10^® 

0.0 

7.2 

U.HO2  +  0  ^  OH  +  O2 

2.00  X  10^® 

0.0 

0.0 

I5.HO2  +  OH  ^  H2O  +  O2 

2.89  X  10^® 

0.0 

-2.1 

Hydrogen  peroxide  reactions 

16^. OH  +  OH  +  M  ^  H2O2  +  M 

kfi  2.30  X  10^® 

-0.9 

-7.1 

koo  7.40  X  10^® 

-0.4 

0.0 

I7.HO2  +  HO2  ^  H2O2  +  O2 

3.02  X  10^2 

0.0 

5.8 

I3.H2O2  +  H  ^  HO2  +  H2 

4.79  X  10^® 

0.0 

33.33 

I9.H2O2  +  H  ^  H2O  +  OH 

1.00  X  10^® 

0.0 

15.0 

2O.H2O2  +  OH  ^  H2O  +  HO2 

7.08  X  10^2 

0.0 

6.0 

2I.H2O2  +  0  ^  HO2  +  OH 

9.63  X  10® 

2.0 

2.0 
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5.2.1  General  Reaction  Rate  Model 


The  reaction  rate  Vj  depends  on  rate  coefficients  and  species  concentra¬ 


tions  42 


(5.23) 


where  kjj  and  kbj  represent  forward  and  backward  rate  coefficients,  which 
are  functions  of  the  temperature  T.  Yk  and  Mk  are  the  mass  fraction  and 
the  molecular  weight  of  the  species  k,  and  and  v'l-  are  reaction  j’s  stoi¬ 
chiometric  coefficients  of  species  k  on  the  reactant  and  product  side. 

The  modeling  of  the  rate  coefficients  kjj  and  k^j  is  an  important  issue 
for  combustion  processes.  The  Arrhenius  law  is  widely  used  to  calculate  the 
forward  rate  coefficient 


(5.24) 


Here  Aj  is  the  pre-exponential  constant,  Uj  the  temperature  exponent,  and 
Ej  the  activation  energy,  all  with  respect  to  reaction  j.  These  constant 
values  are  components  of  a  complete  description  of  the  chemical  scheme  as 
provided  by  the  scientihc  literature,  and  should  be  available  as  input  for 
numerical  computation,  often  in  tabulation  form.  Backward  rate  coefficients 
kbj  are  calculated  from  kfj  and  the  equilibrium  constant  for  reaction  j  by 
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(5.25) 


kbj 


k 


fj 


(P^\ 

\rtJ 


E^=i 


exp  ELi 


ASkj  AH, 


kj 


R  RT 

Here,  po  =  1-0133  x  lO^Pa,  while  AS'fcj  and  Aiffcj  represent  entropy  and 
enthalpy  changes  for  species  j. 

For  reactions  that  involve  a  third  body  effect,  the  reaction  rate  also  de¬ 
pends  on  the  concentrations  of  third  body  components  weighted  by  their 
efficiencies  so  that  the  reaction  rate  coefficient  is  mnltiplied  by  a  factor  cm 


Cm  ^  ^  ViC-iy 
i=l 


(5.26) 


where  rji  is  the  chaperon  efficiency  of  species  i. 

For  some  reactions,  a  more  accnrate  reaction  rate  formnla  is  needed  to 
acconnt  for  the  pressure  and  temperature  dependencies.  In  this  case,  Troe’s 


fitting  formula  43  is  a  popular  approach,  and  the  forward  rate  coefficient  is 
given  by 

k  =  Fkooki,  (5.27) 


where  the  Lindermann-Hinshelwood  form  is 


kr  = 


ko[M]/koo 

1  ko[M]/koo' 


(5.28) 


ko  and  k^o  are  still  computed  by  the  Arrhenius  law,  and  the  corresponding 
parameters  Aq,  A^,  uq,  rioc,  Eq  and  are  still  provided  by  the  mechanism 
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table.  [M]  =  p/RT  is  the  molar  density  of  the  gas  mixture.  A  broadening 
factor  F  to  reduce  the  synthetic  errors  of  the  above  form  of  is  given  by 


logio^ 

N 


_ loglO-^c _ 

1  +  (loglo[^o[^]/^oo]^)^’ 

0.75  -  1.271ogioFc. 


(5.29) 

(5.30) 


The  forms  of  the  pressure  dependent  term  Fc  are  usually  empirical  expres¬ 
sions,  and  can  be  £t  as  a  function  of  temperature  by 


Fe  =  (1  -  a)exp(-T/T”*)  +  ttexp(-T/T*)  +  exp(-T*7T),  (5.31) 


in  which  a,  T*,  T**  and  T***  are  adjustable  parameters  and  should  also  be 
provided  as  part  of  a  mechanisms  table. 

Given  all  the  reactions  in  which  a  certain  species  k  participates,  either  as 
reactant  or  as  product,  we  can  thus  define  the  rate  of  production  of  species 
k  and  the  heat  release  Hvj  of  reaction  j  as 

n  N 

rhk  =  Mk'Yl  i^kjrj,  Hr,j  =  ^  hkMkiykfj,  (5.32) 

j=l  k=l 

where  rj  is  the  reaction  rate  for  reaction  j,  and  hk  is  the  enthalpy  of  species 

k. 
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5.2.2  Reduced  Chemical  Mechanism 


The  chemistry  calculation  contributes  substantially  to  the  cost  of  the  re¬ 
active  flow  simulations.  Analysis  shows  that  when  chemistry  source  terms 
are  calculated  and  provided  as  functions  of  species  concentrations  in  a  system 
of  ordinary  differential  equations  (as  suggested  by  the  hnite  rate  chemistry 
approached  discussed  in  the  previous  chapter),  the  convergence  of  the  chem¬ 
istry  step  is  a  major  computation  bottleneck.  This  is  also  the  case  for  the 
three-dimensional  Scramjet  simulation  discussed  in  Chapter  6.  Decreasing 
both  the  number  of  species  and  the  number  of  reactions  has  the  advantage 
of  signihcantly  decreasing  simulation  complexity. 

We  consider  a  reduction  of  complexity  from  a  global  point  of  view.  In  the 
detailed  mechanism  stated  above,  there  are  eight  reacting  species  involved: 
H2,02,  H2O,  H,0,0H,  HO2  and  H2O2,  and  therefore  eight  corresponding 
differential  equations.  Among  these  species,  H  and  O  are  not  independent, 
they  are  simply  rearranged  among  other  species,  thus  leaving  only  six  inde¬ 
pendent  equations  with  nonzero  chemical  source  terms.  Therefore,  algebraic 
relations  among  hydrogen-oxygen  combustion  system  indicate  that  there  ex¬ 
ists  at  most  six  overall  chemical  steps. 

The  introduction  of  reduced  mechanism  is  proposed  when  it  is  noted 
that  many  elementary  reactions  have  very  little  contribution  to  the  whole 
reaction  process.  There  are  a  number  of  reduced  mechanisms  proposed  in 
the  combustion  literature.  By  analyzing  the  role  and  importance  of  each  of 
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Table  5.3:  Four-step  reduced  mechanism 


(I) 

H  +  02 

OH +  0 

(11) 

0  +  H2  +  M 

H  +  OH +  M 

(III) 

OH +  H2 

H  +  H2O 

(IV) 

H  +  H  +  M 

H2  +  M 

Wi 

Wi  +  We  +  W12  +  Wi5  -  Wi7 
+  Wis  -  W20  +  W21 

Wii 

— 

W2  -  W3  +  We  +  W12  +  wii 
+  W16  +  Wij  -\-  tCig  -|-  W21 

Win 

— 

W3  +  +  Ws  + 

-FW13  -F  Wi5  -1-  Wifi 

Wiv 

— 

W5  +  Wq  +  Wio  —  W12 
+  Wifi  +  Wn  -h  W20 

the  elementary  reactions  of  the  detailed  mechanism,  and  assuming  further 
that  two  intermediate  species  H2O2  and  HO2  are  both  in  steady  state,  a 
four-step  mechanism  has  been  found  to  provide  an  accurate  flame  structure 


for  laminar  diffusion  flames  15  .  The  four  elementary  reactions  chosen  out  of 
the  21  from  the  detailed  mechanism,  along  with  their  algebraic  combinations 


are  given  in  Table  5.3 


The  steady  state  assumptions  are  shown  15  to  provide  accurate  combus¬ 
tion  process,  including  ignition,  extinction  and  heat-release  phenomena  for 
numerical  calculations  for  hydrogen  flames  at  atmospheric  pressure  and  with 
low  boundary  temperatures.  The  coefficients  are  usually  optimized  or  £t  to 
give  correct  key  species  concentration,  maximum  temperature  and  laminar 
premixed  flame  speeds. 
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5.3  Flame  Structures 


5.3.1  Premixed  Flame 


In  a  premixed  flame,  the  fuel  and  oxidizer  are  mixed  perfectly  before  they 
are  injected  into  the  combustion  chamber.  Because  of  its  more  analytically 
explainable  features,  premixed  flames  can  be  regarded  as  a  building  block 
and  a  first  step  towards  the  study  of  more  complicated  flame  structures.  The 


mass  stoichiometric  ratio  defined  as  42 


fYo\  _  ly'oMo 
X^F/st  f'pMf 


and  the  equivalence  ratio  of  the  mixture 


o  In 


st 


(5.33) 


(5.34) 


are  two  central  parameters  commonly  used  to  indicate  the  richness  of  com¬ 
bustion  regimes,  with  0  >  1  for  rich  combustion,  when  the  fuel  is  in  excess 
and  (f)  <1  for  lean  combustion,  when  the  oxidizer  is  in  excess.  Here,  Yp  and 
Yq  stand  for  the  mass  fractions  of  fuel  and  oxidizer,  respectively. 

The  determination  of  flame  thickness  is  crucial  to  the  numerical  simula¬ 
tion  of  reacting  flows.  The  flame  structure  should  be  resolved  with  enough 
cells  located  in  the  flame  front.  There  are  several  ways  to  define  the  width 
of  a  flame  <5°.  For  numerical  purposes,  a  practical  formula  is  proposed  from 
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the  temperature  profile  42 


5? 


T  —  T  • 

max  mm 


max 


(lii) 


(5.35) 


5.3.2  Diffusion  Flame 

In  a  diffusion  flame,  the  fuel  and  oxidizer  are  not  premixed.  They  are 
added  into  a  combustion  chamber  separately  through  different  inlets.  In  dis¬ 
tinction  to  a  premixed  flame,  the  mixing  is  characterized  by  a  global  equiv¬ 
alence  ratio 


where  nip  and  rrio  are  the  flow  rate  of  fuel  and  of  oxidizer  respectively. 

Compared  to  a  premixed  flame,  diffusion  plays  a  much  more  important 
role  here.  The  chemical  reactions  take  place  primarily  near  the  area  where 
fuel  and  oxidizer  are  not  far  from  the  stoichiometric  ratio.  Burning  is  far 
less  efficient  than  for  a  premixed  flame,  considering  the  lack  of  either  fuel 
or  oxidizer  on  the  far  end  of  both  sides.  Furthermore,  the  flame  speed  and 
the  thickness  of  the  diffusion  flame  are  not  as  well-dehned  as  for  premixed 
flames,  the  former  due  to  the  inability  of  the  flame  to  propagate  towards  the 
fuel,  and  the  latter  due  to  the  slow  yet  steady  growth  of  the  reaction  zone 
as  the  time  progresses.  However,  when  a  hxed  time  frame  is  considered  (a 
flow  through  from  the  fuel  inlet  to  the  outlet  of  the  Scramjet  combustor  for 
example),  the  range  of  the  reaction  zone  can  still  be  dehned.  In  fact,  the 
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full  width  at  half  maximum  (FWHM)  of  the  temperature  prohle  is  widely 
considered  to  be  a  characteristic  flame  thickness  of  diffusion  flames.  This 
dehnition  is  used  in  the  present  context  to  study  the  diffusion  flame  and  its 
mesh  requirements. 


5.3.3  Thickened  Flame  Model 

The  thickened  flame  model  is  a  method  developed  to  thicken  the  flame 
front  artihcially,  to  allow  numerical  solution  of  the  flame  front  on  a  coarse 
grid.  The  level  of  increased  thickness  is  achieved  by  decreasing  the  pre¬ 
exponential  constant  A  of  the  previously  discussed  Arrhenius  form  by  a  fac¬ 
tor  of  TF,  and  increasing  by  the  same  factor  the  level  of  the  species  and 
heat  diffusion  coefficients.  Theories  of  laminar  flame  structure  show  that  the 
expression  of  flame  speed  and  flame  thickness  scale  like 


S  r  oc 


DjhRr, 


sr 


oc 


Dth 
.0  ’ 


(5.37) 


where  and  Rr  are  the  diffusion  coefficient  and  mean  reaction  rate,  respec¬ 
tively.  Therefore,  by  replacing  Dth  by  TFDth  and  A  by  A/F,  the  reaction 
rate  Rr  is  replaced  by  Rr/TF,  the  thickened  flame  propagates  at  the  same 
laminar  flame  speed  sl,  but  with  a  flame  width  thickened  to  TFS^.  Choosing 
a  sufficiently  large  TF  (between  10  to  100  in  most  cases)  allows  the  flame  to 
be  resolved  on  an  LES  grid  P. 


It  is  shown  in  Section  5.5.3  that  thickened  flame  model  tends  to  produce 
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increased  errors  when  modeling  diffusion  ffames.  In  fact,  the  thickened  flame 
model  was  originally  derived  for  premixed  ffames.  To  adapt  thickened  flame 
to  the  case  of  diffusion  ffames  and  more  complicated  flame  regimes,  a  new 
model  to  replace  the  constant  thickening  factor  with  a  local  variable  that 
depends  on  space  and  time  was  proposed  by  J.  P.  Legier  et.al  in  [^.  As 
the  name  suggests,  the  dynamically  thickened  flame  (DTP)  model  adjusts 
the  thickening  factor  dynamically  across  the  whole  region.  The  idea  is  to 
maintain  the  specified  thickening  level  inside  the  combustion  zone,  while 
suppressing  it  away  from  the  flame  front.  In  this  way,  mixing  is  correctly 
represented  outside  the  reaction  zone  and  the  flame  itself  still  maintains 
enough  thickness  to  be  resolved. 

In  its  original  form,  an  indicator  used  to  determine  the  range  of  combus¬ 
tion  is  defined  by 

f!  =  y;'yyexp  (-rP I ,  (5.38) 

where  is  the  activation  temperature  and  T  (>  1)  is  a  specified  parameter 
that  further  decreases  the  activation  temperature.  The  sensor  is  then  applied 
to  obtain  the  variable  thickening  factor  in 

TF  =  1 -h  (TFmax  -  l)tanh  (5.39) 

Here,  Umax  is  the  maximum  possible  value  of  and  /3  is  another  parameter 
that  controls  the  thickness  of  the  transition  layer  between  thickened  and  non- 
thickened  zones  [^.  In  the  present  analysis,  we  replace  the  indicator  with 
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the  local  reaction  rate  which  serves  the  same  purpose.  Simulation  results 
show  that  the  major  flame  properties  (flame  thickness,  mean  concentration 
of  OH  and  HO2,  ignition  delay,  etc.)  stay  at  the  same  level  as  when  applying 
the  original  thickened  flame  model  using  the  same  thickening  values.  This 
implies  that  DTF  retains  the  effect  of  the  TF  model  within  the  flame  region, 
with  the  additional  benefit  that  outside  of  the  region,  diffusion  rates  are 
not  affected  by  the  model.  Moreover,  the  parameter  {3  provides  increased 
flexibility  in  adjusting  the  thickening  at  the  edge  of  the  flame  front. 


5.4  Combustion  Models 


5.4.1  The  Flamelet  Model 

To  derive  the  model,  we  begin  with  the  continuity  equation 


dp 


+ 


dpuj 

dxi 


0, 


(5.40) 


where  p  is  the  density,  Ui  is  the  velocity  component,  and  t  is  time.  When 
dealing  with  a  mixture  of  several  species,  a  continuity  equation  is  required 
for  each  of  them.  In  the  species  equation,  an  extra  term  appears  due  to  the 
rate  of  production  of  that  species  by  chemical  reactions.  This  quantity  is 
denoted  by  m*.  Therefore,  the  continuity  equation  for  the  partial  density  pi 
of  species  i  can  be  written  as 


dpYk 

dt 


+ 


A 

dxi 


{p{ui  +  Vk,i)Yk)  =  rhfc. 


(5.41) 


46 


where  pi  and  refer  to  partial  density  and  diffusion  velocity. 
After  applying  Pick’s  law,  the  above  can  be  rewritten  as 


dpYk  dpYkUj 
dt  dxi 


A 

dxi 


=  rfik- 


(5.42) 


The  key  idea  of  the  flamelet  approach  is  a  change  of  variables  within 
the  diffusion  flame  sheets.  The  underlying  assumption  is  that  the  structure 
of  the  flame  depends  on  the  mixture  fraction  Z  and  time  t  only,  so  that 
temperature  and  species  mass  fractions  can  be  expressed  as  T  =  T{Z,t)  and 
Yk  =  Yk{Z,t).  Therefore,  when  a  variable  change  is  exercised  on  species 
equations  from  {xi,X2,X3,t)  to  {Z,yi,y2,t),  the  gradient  terms  with  respect 
to  the  directions  of  yi  and  y2  can  be  neglected  when  comparing  to  the  terms 
normal  to  the  flame.  Here  yi  and  y2  are  spatial  variables  orthogonal  to  each 
other  in  the  plane  normal  to  direction  defined  by  Z.  The  assumption  implies 
that  X  is  represented  uniquely  by  the  mixture  fraction  coordinate  Z,  and 
that  the  reaction  zone  is  smaller  than  the  small  scales  of  the  turbulence. 
The  flamelet  derivation  is  restricted  to  a  small  region  of  locally  laminar  flow 
around  the  reaction  zone. 

After  converting  to  local  coordinates  {Z,  t)  and  rearranging  terms 


P 


dYk 

dt 


+  Yk 


+ 


dYk 

dZ 

pD^ 


(  dZ  dZ  d 

At  _  • 

dxi  dxi )  dx"^ 


(5.43) 
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The  factors  within  the  parentheses  of  the  second  and  third  terms  on  the  left 


of  Equation  |5.43|  vanish  due  to  conservation  of  mass  and  species  equation, 
and  the  equation  becomes 


dY. 


P- 


dt 


=  mk  +  pD 


dZ  dZ  \  d'^Y, 


dxi  dxi  J  dZ^ 


1  d^Yk  . 


(5.44) 


In  this  way,  the  flamelet  model  effectively  reduces  the  system  of  species  equa- 

,  'dZ  dZ 

tions  to  one  dimension,  the  so-called  Z-space.  The  variable  y  =  2D 


dxi  dx 

is  called  the  scalar  dissipation  rate.  It  is  a  key  variable  in  the  flamelet 
equation;  it  controls  mixing  and  measures  molecular  fluxes  between  differ¬ 
ent  species  in  directions  normal  to  the  flame  front.  Once  y  is  specihed,  the 
combustion  problem  is  solved  and  the  flame  structure  provided.  In  other 
words,  the  hydrodynamics  has  been  removed  from  the  chemistry  in  Equation 


5.44  which  is  a  pure  reaction  diffusion  equation.  The  steady  flamelet  model 
assumes  further  that  the  flame  structure  is  in  steady  state,  and  the  flamelet 
equation  can  be  simplihed  to 


1  d^Yk 

2™^  +  ™  =  0- 


(5.45) 


One  approach  to  simulating  the  Scramjet  flow  structure  is  to  use  the 
flamelet  model  described  above.  The  model  is  computationally  efficient  in 
that  the  ffamelet  solutions  can  be  precomputed  and  tabulated  for  further  use. 
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given  specific  boundary  conditions.  These  tables  describe  scalar  variables  as 
functions  of  Z  and  y.  When  the  flow  equations  are  solved  in  main  simulation, 
a  table  lookup  is  then  conducted  for  each  cell. 


5.4.2  Finite  Rate  Chemistry 


As  Section  |5.4.1|  suggests,  one  signihcant  assumption  inherent  in  the 
flamelet  model  is  that  the  flame  is  considered  to  be  a  locally  one-dimensional 
structure,  and  quantities  that  account  for  gradients  parallel  to  the  flame  front 
are  therefore  ignored.  A  numerical  scheme  based  on  this  assumption  may  not 
accurately  model  flame  conhgurations  in  many  practical  problems,  especially 
when  turbulence  plays  a  major  role  in  the  flame  development. 

Finite  rate  chemistry  is  a  method  developed  here  to  characterize  turbulent 
combustion  processes  as  an  alternative  to  the  flamelet  model.  It  is  based  on 
the  notion  that  turbulence  occurs  on  many  different  levels  of  length  scales, 
and  that  the  chemistry  time  and  length  scales  he  somewhere  between  the 
smallest  turbulence  length  scale  (the  Kolmogorov  scale)  and  the  larger  end 
of  it.  Therefore,  it  is  possible  that  combustion  process  can  be  resolved  on  a 
larger  scale  than  the  levels  needed  to  resolve  turbulence.  As  an  advantage,  the 
model  does  not  rely  on  the  laminar  flame  structure  assumption.  Instead,  the 
combustion  processes  and  rate  of  chemical  reactions  are  calculated  directly 
with  data  provided  by  chemical-kinetic  mechanisms.  A  system  of  ordinary 


differential  equations  as  shown  in  (5.46 )  based  on  the  reacting  species  are  then 


solved  at  each  time  step  to  account  for  the  change  in  species  concentration 
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due  to  the  combustion  process. 


m 

dt 


=  mk  = 


Mk'Y^h’kjTj.  {k  =  l,...,N) 


(5.46) 


i=i 


5.5  Finite  Rate  Chemistry  Modeling  in  One  Dimen¬ 


sion 


In  this  section,  the  hnite  rate  chemistry  model  discussed  in  Section  |5.4.2 
is  investigated  in  a  one-dimensional  framework.  The  reduction  of  chemical- 
kinetic  mechanism  and  thickening  factors  are  both  related  to  the  quality  of 
simulation  results.  The  purpose  of  the  investigation  is  to  analyze  the  accuracy 
of  these  approaches  and  to  provide  necessary  guidance  to  the  mesh  resolution 
needed  to  perform  a  full-domain  numerical  simulation  of  a  three-dimensional 
Scramjet  turbulent  combustion  problem. 


5.5.1  Simulation  Setup 

To  simulate  one- dimensional  combustion,  various  formulations  can  be  de¬ 
veloped  to  study  different  types  of  flames  34,^^.  In  this  context,  initial 
conditions  are  set  up  consistent  with  the  three-dimensional  Scramjet  flow 
conhgurations.  Densities,  temperatures,  velocities  and  species  partial  densi- 
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ties  are  set  to  be 


X  =  0  :  pL  =  pH^  =  1.241  X  10■^  Tl  =  SOO/sT,  vl  =  0,  (5.47) 

x  =  L:  pp.  =  po^  =  1.241  X  10■^  Tr  =  150077,  vr  =  0.  (5.48) 


The  partial  densities  are  the  conhguration  of  a  diffusion  flame.  To  study  a 

premixed  type  of  flame,  a  stoichiometric  mixture  of  hydrogen  and  oxygen  is 

1  s 

considered,  so  that  pho  = - P  and  po,  = - p  on  both  ends. 

n  2  1  +  1  +  s 

A  solution  for  this  boundary  value  problem  can  be  found  by  discretizing 
the  spatial  domain  with  appropriate  numerical  schemes.  The  WENO  scheme 
as  described  in  Section  PI  combined  with  finite  rate  chemistry  discussed 
in  Section  |5.4.2  are  applied  here  to  study  the  flame  structure,  which  is  a 
simplified  version  of  the  combustion  process  taking  place  inside  the  Scramjet 
combustor.  Both  the  detailed  chemical-kinetic  mechanism  in  Table  15.  II  and 


the  reduced  version  in  Table  5.3  are  employed.  In  addition,  it  is  assumed 
that  all  binary  diffusion  coefficients  are  constant  and  equal:  D  =  Dk  =  10“^. 
Viscosity  and  thermal  conductivity  are  determined  according  to  Section  5.1[ 
The  simulation  consists  of  five  parts:  the  chemistry  data  input  (chemistry- 
kinetic  mechanism  table  and  thermodynamic  table),  the  reaction  table  parser 
and  species  and  reaction  data  structures,  the  general  numerical  scheme,  the 
equation  of  state  library,  and  the  hnite  rate  chemistry  ordinary  differential 
equation  solver. 
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5.5.2  The  One  Dimensional  Flame  Structure 


A  typical  simulation  result  of  the  one-dimensional  flame  structures  for 


premixed  and  diffusion  flames  is  illustrated  in  Figure  5.1 


The  structure  of  a  premixed  flame  can  be  roughly  divided  into  four  parts: 
the  unburned  mixture,  the  preheat  zone,  the  reaction  zone  and  the  burned 
mixture.  In  unburned  and  burned  zones,  temperature  and  concentration  of 
fuel/oxidizer  are  almost  constant.  In  the  preheat  zone,  temperature  starts 
to  grow  due  to  the  amount  of  net  heat  release  contributed  by  exothermic 
reactions.  It  is  the  region  where  the  unburned  gas  mixture  is  heated  up  by 
conduction.  The  reaction  zone  follows  the  preheat  zone,  and  is  the  region 
where  major  reactions  take  place.  The  gas  mixture  coming  out  of  the  re¬ 
action  zone  then  enters  the  burned  gas  zone  where  temperature  reaches  a 


high  level  and  stays  constant  again.  Also  shown  in  Figure  5T  are  the  mass 
concentration  prohles  of  the  fuel  ^^2,  the  oxidizer  O2,  and  OH,  a  product 
species  of  major  importance  in  hydrogen-oxygen  combustion. 

The  grid  resolution  needed  to  resolve  a  hydrogen-oxygen  premixed  flame 
is  generally  considered  impractical  for  most  simulation  purposes,  with  a  flame 


width  roughly  in  the  order  of  between  10“^  cm  to  10“^  cm  [^[8,27  ,  meaning 
that  the  mesh  size  for  a  well-resolved  combustion  scheme  should  be  at  least  in 
the  order  of  10“^  cm  (to  guarantee  about  10  cells  inside  of  the  flame  front). 
This  range  of  flame  widths  is  verihed  by  the  one-dimensional  simulations 
performed  here  at  a  mesh  size  of  0.001cm  to  about  0.05  cm,  using  the  flame 
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Figure  5.1:  The  flame  structure  of  premixed  (upper)  and  diffusion  (lower) 
flames.  The  variables  included  are  temperature  and  the  mass  fractions  of  the 
oxidizer  (O2),  the  fuel  {H2)  and  one  indicator  product  species  OH. 
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width  formula  given  by  Equation  |5.35 

As  Figure  0  shows,  a  diffusion  flame  exhibits  a  very  different  flame 
structure.  The  combustion  is  dominated  by  the  rate  of  mass  diffusion  between 
the  fuel  {H2)  and  the  oxidizer  (O2).  Generally  speaking,  the  burning  is  far 
less  efficient,  as  the  close  to  unity  mass  fractions  of  fuel  and  oxidizer  away 
from  the  flame  zone  indicate.  On  the  other  hand,  it  can  be  seen  that  the 
mass  concentration  of  species  OH  is  a  good  indicator  of  the  temperature 
level  in  both  premixed  and  diffusion  flames. 

Due  to  the  effect  of  mass  diffusion,  the  width  of  a  diffusion  flame  is  gen¬ 
erally  much  thicker  than  that  of  a  premixed  flame.  For  this  reason,  modeling 
a  diffusion  flame  is  less  restrictive  in  terms  of  grid  resolution  requirement. 
However,  when  the  thickened  flame  model  is  applied  to  a  flame  with  a  more 
complex  regime,  as  is  the  case  of  a  flame  inside  a  combustion  chamber,  the 
flame  prohle  and  the  combustion  process  may  be  more  susceptible  to  a  change 
caused  by  the  thickening  factor,  and  therefore  the  flame  structure  may  not 
be  correctly  predicted. 

The  effectiveness  of  the  thickened  flame  model  and  its  influence  on  the 
key  properties  of  the  flame  structure  are  discussed  below. 

5.5.3  Thickened  Flame  and  Grid  Resolution 

In  this  section,  we  investigate  the  effect  of  the  thickened  flame  model 
and  the  level  of  reduction  in  chemical  mechanism  on  the  quality  of  the  one¬ 
dimensional  flame  structure  by  performing  a  list  of  tests  that  are  aimed  at 
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capturing  the  major  flame  characteristics.  A  combination  of  model  parame¬ 
ters  that  give  the  correct  flame  predictions  are  to  be  derived. 

As  discussed  in  Section  5.5.2[  premixed  flames  are  generally  hard  to  re¬ 
solve.  The  thickened  flame  model  can  therefore  be  employed  to  reduce  the 
mesh  requirements  needed  to  resolve  a  premixed  flame.  A  flame  is  generally 
considered  to  be  resolved  if  certain  indicator  parameters  of  the  flame  are 
preserved  and  shown  to  be  insensitive  to  changes  in  discretization  size  (and 
thus  are  converging  to  the  true  value).  In  this  section,  three  such  indicator 
parameters  are  investigated,  namely,  the  laminar  flame  speed,  the  peak  mass 
fraction  of  OH  (a  major  product  of  general  interest)  and  the  peak  mass  frac¬ 
tion  of  HO2  (a  very  active  intermediate  radical  species).  The  last  quantity 
is  generally  considered  to  be  difficult  to  achieve  convergence  for  because  of 
its  unstable  nature  and  due  to  the  fact  that  it  is  produced  in  considerable 
amount  only  in  a  very  thin  layer  of  reaction  zone,  which  is  coincident  with 
the  location  of  the  flame  front. 

On  the  other  hand,  a  diffusion  flame  maintains  a  very  different  structure, 
and  the  flame  speed  is  not  well-defined.  In  addition  to  this,  species  mass 
productions  are  not  conserved  quantities  across  different  levels  of  mesh  size. 
Due  to  the  extent  of  mass  diffusion,  finer  grids  generally  contribute  to  more 
burning,  and  therefore  higher  proportions  of  product  species.  What  turns 
out  to  be  an  important  indicator  for  a  diffusion  flame  is  the  delay  in  ignition 
time.  The  ignition  time  of  a  flame  regime  can  be  roughly  divided  into  the 
time  delay  due  to  mixing  and  the  delay  for  the  activation  energy  to  be  reached 
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and  the  chemical  reactions  to  occur.  While  the  hrst  type  of  delay  is  almost 
negligible  for  a  premixed  flame,  it  could  be  quite  signihcant  for  a  diffusion 
flame.  Moreover,  the  combined  effects  of  thickened  flame  model  (increases 
diffusion  while  decreasing  reaction  rates)  tend  to  further  delay  the  ignition, 
and  the  amount  of  heat  release  produced  by  combustion.  Therefore,  the 
ignition  time  is  instead  studied  as  a  criteria  in  the  case  of  diffusion  flames. 


5. 5. 3.1  Premixed  Flame  The  laminar  flame  speed  is  computed  using 
a  number  of  grid  sizes  ranging  from  Ax  =  0.001cm  to  Ax  =  0.09  cm.  The 
thickened  flame  model  is  applied  on  three  levels:  TF  =  1  (i.e.  no  thicken¬ 
ing),  TF  =  5  and  TF  =  10.  Both  full  and  reduced  chemical  mechanisms 
are  considered.  The  objective  is  to  hnd  the  coarsest  mesh  size  in  each  com¬ 
bination  of  TF  and  chemistry  model  that  predicts  the  flame  speed  to  an 
acceptable  degree.  The  flame  speed  and  species  mass  fractions  calculated 
under  the  hnest  mesh  size  Ax  =  0.001cm  are  considered  as  benchmark  val¬ 
ues.  When  different  levels  of  thickening  factor  are  imposed,  the  flame  speed  at 
this  hnest  grid  stays  at  around  76cm/ms,  the  mass  fraction  Yho2  ~  0.325% 
and  Yqh  ~  14.46%,  for  both  full  and  reduced  mechanisms.  The  thickness  of 
hame  is  evaluated  to  be  about  0.04  cm,  a  level  clearly  not  resolvable  under 
the  current  Scramjet  mesh  size  of  0.05  cm.  Using  a  maximum  allowed  10% 
error  in  the  hame  speed  and  the  mass  fraction  of  species  as  the  acceptance 


level,  results  are  shown  in  Figure  5.2 


The  result  indicates  a  few  points. 
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Full  Mechanism  TF  =  1 


Reduced  Mechanism  TF  =  1 


Mesh  Size  (cm)  Y_FI02  Mesh  Size  (cm)  Y_OH 

Full  Mechanism  TF  =  5  Reduced  Mechanism  TF  =  5 


Mesh  Size  (cm)  Y_H02 


Full  Mechanism  TF  =  10 


Mesh  Size  (cm)  Y_HO 
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Figure  5.2:  Simulated  premixed  flame  speed  and  species  concentration  {HO2 
for  full  chemistry  and  OH  for  reduced  chemistry)  against  a  selection  of  mesh 
sizes.  Full  chemistry  (left  column)  and  reduced  chemistry  (right  column) 
with  no  thickening  TF  =  1  (row  1),  TF  =  5  (row  2),  and  TF  =  10  (row  3). 
On  these  figures,  the  dotted  horizontal  lines  indicate  the  correct  chemistry 
while  the  solid  horizontal  line  indicates  the  acceptance  level. 


57 


1.  When  thickened  flame  is  not  applied  {TF  =  1),  modeling  combustion 
accurately  requires  a  very  high  level  of  mesh  refinement,  in  the  order  of 
10“^  cm  for  full  mechanism  and  around  10“^  cm  in  the  case  of  reduced 
mechanism; 

2.  The  thickened  flame  model  effectively  expands  the  width  of  the  flame 
front,  making  it  possible  to  resolve  the  flame  width  on  a  much  coarser 
grid.  Moreover,  at  least  for  small  thickening  ranges,  the  allowed  coars¬ 
est  mesh  size  is  roughly  proportional  to  the  level  of  thickening; 

3.  A  reduced  chemical  mechanism  allows  a  much  coarser  grid  size  when 
it  comes  to  modeling  the  flame  speed.  The  reason  is  that  reduced 
mechanisms  are  generally  fit  to  the  detailed  mechanism  in  terms  of  the 
laminar  flame  speed.  Fewer  reactions  with  comparable  reaction  rates 
and  the  absence  of  unstable  intermediate  species  also  contribute  to  a 
more  stable  flame  speed  in  contrast  to  the  detailed  mechanism; 

4.  The  intermediate  species  HO2  is  assumed  to  be  in  a  steady  state  when 
applying  the  reduced  mechanism.  Therefore  in  this  case,  OH  is  used 
as  a  replacement  in  the  reduced  mechanism  plots.  It  becomes  obvious 
from  the  plots  that  OH  is  a  more  stable  product  species  comparing  to 
HO2,  and  it  may  be  less  indicative  of  the  correct  level  of  mesh  size  that 
can  be  used. 


The  mesh  requirement  results  are  summarized  in  Table  5.4 
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Table  5.4:  Premixed  flame  mesh  requirements  with  respect  to  full  and  re¬ 
duced  chemistry  mechanism  and  thickening  factors  TF  =  1,  TF  =  5  and 
TF  =  10.  The  mesh  resolutions  are  based  on  a  requirement  for  a  maximum 
10%  error  in  flame  speed  and  species  mass  fractions. 


Full  mechanism 

Reduced  mechanism 

TF  =  1 

TF  =  5 

TF  =  10 

TF  =  1 

TF  =  5 

TF  =  10 

0.004  cm 

0.025  cm 

0.05  cm 

0.01  cm 

>  0.05  cm 

0.07  cm 

One  disadvantage  of  the  thickened  flame  model  is  that  the  artificial  thick¬ 
ening  may  alter  other  flame  properties.  Due  to  thickening,  the  interaction 
between  the  flame  and  the  turbulent  flow  structure  may  be  affected,  espe¬ 
cially  flow  structure  on  small  scales.  To  reduce  this  effect,  only  the  smallest 
possible  thickening  level  should  be  used.  The  choice  of  thickening  level  is 
therefore  not  arbitrary.  It  is  suggested  in  numerous  literature  (for  example 
in  that  a  thickening  factor  can  take  on  values  as  large  as  100,  however, 
the  study  below  shows  that  such  extremes  are  both  undesirable  (for  instance. 


causing  the  delay  in  ignition  time  as  described  in  Section  5. 5. 3. 2)  and  unneo 


essary  for  the  present  problem.  Actually  a  limited  thickening  is  sufficient  to 
resolve  the  flame  front  at  least  in  the  flame  configuration  studied  here. 

The  desired  value  of  thickening  factor  to  be  applied  has  been  considered 
by  others.  For  instance,  proposes  that  a  suitable  value  TF^ax  could  be 


_ ....  (nAjnesh  .  , 

^  ^max  niSiX  ^  ,  1  J  , 


(5.49) 


where  A^esh  is  the  computational  mesh  size,  6i  is  the  flame  width,  and  n  rep¬ 
resents  the  number  of  cells  inside  of  the  flame  front  to  resolve  the  thickness. 
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n  is  usually  chosen  in  a  empirical  way,  and  5  is  considered  an  appropriate 
value.  Instead,  here  we  choose  to  determine  the  thickening  level  following  the 
principle  that  the  artihcially  thickened  flames  should  still  be  able  to  maintain 


the  macro  flame  structure  and  combustion  properties.  Figure  5.3  again  gives 
the  simulated  results  of  the  flame  speed  and  species  concentration  of  HO2  for 
full  chemical  mechanism  and  OH  for  reduced  chemical  mechanism.  A  range 
of  thickening  factors  from  TF  =  ltoTF  =  20  are  applied  to  determine 
the  smallest  acceptable  value  of  TF,  given  a  range  of  mesh  sizes  between 
0.01cm  and  0.05  cm.  Figure  |5]3  shows  the  cases  for  Ax  =  0.05  cm,  which  is 
consistent  with  the  current  Scramjet  mesh  resolution,  and  Ax  =  0.025  cm, 
which  is  a  further  rehned  mesh  size  that  can  be  used. 

The  result  indicates  that  the  coarser  grid  generally  requires  more  thick¬ 
ening,  as  expected.  The  number  of  cells  in  the  flame  front  determined  in 
this  way  are  6  for  the  grid  level  of  Ax  =  0.05  cm  and  8  for  Ax  =  0.025  cm, 


more  demanding  than  given  in  (5.49).  Also,  the  two  parameters  considered 
both  tend  to  be  less  sensitive  to  an  increase  in  thickening  level.  This  does 
not  mean  that  the  flow  features  are  correctly  captured.  To  maintain  a  flame 
structure  as  close  to  the  real  flame  as  possible,  only  the  smallest  TF  value 
should  be  applied.  In  addition,  the  plots  on  the  reduced  mechanism  side 
indicate  a  constant  over  prediction  of  the  species  concentration  of  OH .  This 
is  an  error  inherent  in  the  development  of  reduced  mechanism,  a  topic  to  be 
discussed  further  in  Section  15.5.3.31 

To  summarize,  a  premixed  flame  needs  to  be  thickened  in  most  applica- 
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Reduced  Mechanism  (mesh  size  =  0.05  cm) 
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Figure  5.3:  Simulated  premixed  flame  speed  and  species  concentration  {HO2 
for  full  chemistry  and  OH  for  reduced  chemistry)  against  a  selection  of  thick¬ 
ening  factors.  Full  chemistry  (left  column)  and  reduced  chemistry  (right  col¬ 
umn)  with  Ax  =  0.05  cm  (row  1)  and  Ax  =  0.025  cm  (row  2)  are  displayed. 
The  meanings  of  dotted  and  solid  horizontal  lines  have  the  same  meaning  as 


is  described  in  Figure  5.2 
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Table  5.5:  The  minimum  level  of  thickening  factor  needed  to  resolve  a  pre¬ 
mixed  flame  for  grid  resolution  ranging  from  Ax  =  0.01  cm  to  Ax  =  0.05  cm. 
Both  the  full  and  the  reduced  chemical  mechanisms  are  considered. 


Full  chemical  mechanism 

0.01  cm 

0.02  cm 

0.025  cm 

0.03  cm 

0.04  cm 

0.05  cm 

TF  ■  —  2 

mm 

TF  —  2 

TF  —  ^ 

mm  ^ 

T  F  ■  =  6 

T  F  ■  =  6 

T  F  '  =8 

mm  ^ 

Reduced  chemical  mechanism 

0.01  cm 

0.02  cm 

0.025  cm 

0.03  cm 

0.04  cm 

0.05  cm 

TF  ■  —  1 

^  mm 

TF  ■  —  1 

TF  —  2 

mm  ^ 

TF  ■  —  2 

^  mm  ^ 

TF  ■  —  2 

mm  ^ 

TF  —  4 

mm  ^ 

tions,  when  mesh  requirements  are  not  met.  On  the  other  hand,  thickening 
should  not  be  introduced  for  more  than  needed.  The  minimum  level  of  thick¬ 
ening  factor  for  both  full  and  reduced  chemical  mechanisms  are  studied,  and 


the  concluding  results  for  premixed  flames  are  summarized  in  Table  5.5 


5. 5. 3. 2  Diffusion  Flame  Based  on  the  results  obtained  by  Section 


5.5.3.1 


this  section  goes  on  to  examine  the  mesh  requirements  for  correctly  simu¬ 


lating  a  one-dimensional  diffusion  flame.  As  explained  in  Section  5.3.2,  the 
thickness  of  a  diffusion  flame  is  dehned  in  a  different  way.  Because  of  its 
slowly  developing  nature,  a  steady  flame  front  in  the  normal  sense  does  not 


exist,  unless  the  flame  itself  is  stretched  and  certain  conditions  met  42  .  In 


this  context,  we  restrict  the  combustion  time  period  to  T  =  0.12  ms,  which 
is  approximately  the  time  needed  for  the  flow  of  fuel  {H2)  to  travel  from  the 
jet  nozzle  to  the  end  of  the  Scramjet  combustion  chamber,  and  therefore  the 
maximum  time  allowed  for  a  diffusion  type  of  flame  to  develop  inside  the 
chamber.  The  thickened  flame  model  is  applied  in  conjunction  with  both  the 
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full  and  reduced  chemical  mechanisms.  The  flame  width  dehned  by  FWHM 


described  in  Section  5.3.2  turns  out  to  be  much  thicker  than  that  of  premixed 
flame  under  the  same  level  of  mesh  resolution.  The  hnest  mesh  in  this  study 
(Ax  =  0.001cm)  evaluates  the  diffusion  flame  thickness  to  be  about  0.2  cm 
at  the  end  of  the  time  period. 


It  is  suggested  in  Section  |5.3.2|  that  a  diffusion  flame  is  usually  located 
where  the  mixture  of  fuel  and  oxidizer  is  close  to  a  stoichiometric  ratio,  and 
the  extent  of  mixing  is  an  important  factor  in  determining  the  efficiency  of 
diffusion  flame  combustion.  Therefore,  artihcial  thickening  may  affect  the 
diffusion  flame  structure  more  negatively  comparing  to  a  premixed  flame,  in 
view  of  the  fact  that  the  interactions  between  mass  diffusion  and  chemical 
reactions  are  highly  nonlinear.  One  of  the  major  signihcant  influences  is  on 
the  ignition  time  of  the  combustion  process.  A  delay  of  ignition  to  a  certain 
degree  could  affect  the  flame  prohle  as  a  whole  and  the  extent  of  combustion, 
especially  when  the  time  period  for  combustion  to  take  place  and  for  the  flame 
to  develop  is  limited,  as  is  the  case  here. 

In  the  thickening  analysis  of  the  premixed  flames,  we  restrict  the  range  of 
T F  that  could  be  applied  to  a  series  of  mesh  sizes  ranging  from  Ax  =  0.01  cm 
to  Ax  =  0.05  cm  in  order  to  obtain  simulation  results  within  an  acceptance 
level  of  10%.  Based  on  the  result  of  minimum  TF  values  in  each  cases  as 


summarized  in  Table  |5.5[  we  further  determine  in  each  case  whether  the 
minimum  level  of  TF,  when  imposed  on  a  diffusion  flame,  causes  too  much 
ignition  delay. 
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Full  Mechanism  Ax  =  0.05  cm 


Reduced  Mechanism  Ax  =  0.05  cm 


Reduced  Mechanism  Ax  =  0.025  cm 


Figure  5.4:  The  evolution  of  diffusion  flame  heat  release  over  the  time  period 
T  =  0.12  ms.  Thickening  factors  derived  from  the  premixed  flame  analysis 
are  applied  in  each  case  (left  column  for  full  mechanism,  and  right  column 
for  reduced  mechanism). 


Figure  |5~4|  illustrates  the  evolution  of  generated  heat  release  over  the  time 
period  T  =  0.12  ms,  for  two  levels  of  grid  size  Ax  =  0.05  cm  to  Ax  =  0.025  cm 
and  two  levels  of  chemical  mechanism. 

It  can  be  seen  from  the  hgure  that  applying  the  thickened  flame  model 
generally  postpone  the  flame  ignition  time.  Between  the  two  grid  levels 
investigated  here,  the  coarse  grid  (Ax  =  0.05  cm)  shows  a  larger  discrepancy 
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between  no  thickening  {TF  =  1)  and  TF  starts  at  8  (the  minimum  value 
determined  by  premixed  flame).  The  ignition  is  not  started  until  a  little 
more  than  half  of  allowed  burning  time  has  passed.  Due  to  the  delay,  the 
total  heat  release  at  the  end  of  simulation  time  is  reduced  by  about  28%. 
With  the  reduced  mechanism,  where  TF  is  justified  to  start  from  a  smaller 
value,  the  delay  is  about  the  same  extent,  though  after  the  point  of  ignition 
is  passed,  the  rate  of  heat  release  increases  rapidly,  with  an  end  heat  release 
error  stopping  at  around  8%. 

With  one  more  level  of  mesh  refinement,  the  difference  introduced  by 
thickening  factor  narrows  by  a  lot,  partly  due  to  the  observation  that  a 
smaller  TF  could  be  applied  with  more  mesh  refinement.  Considering  the 
minimum  TF  levels  for  instance,  the  ignition  delay  amounts  to  29%  and  15%, 
the  reduction  in  total  heat  release  amounts  to  9.6%  and  —0.01%  for  the  full 
and  reduced  chemical  mechanisms,  respectively. 


The  dynamically  thickened  flame  model  introduced  in  Section  5.3.3  shows 
no  large  effect  on  flame  features  in  the  one-dimensional  diffusion  flame  con¬ 
figurations.  An  analysis  of  the  result  shows  that  there  is  very  little  chemical 
reaction  taking  place  outside  of  the  reaction  zone  even  without  applying  the 
dynamic  model.  This  indicates  that  increased  diffusion  by  thickened  flame 
model  has  very  little  effect  on  the  combustion  profiles  there.  However,  com¬ 
bustion  processes  in  the  three-dimensional  framework  is  more  complicated, 
with  much  more  interaction  between  different  types  of  flame  structures.  The 
dynamically  thickened  flame  model  could  be  more  preferable  due  to  its  ability 
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to  maintain  the  correct  mixing  profile. 

The  production  of  OH  and  HO2  species  discussed  earlier  turns  out  to  be 
more  sensitive  to  a  change  in  thickening  factor  in  a  diffusion  flame  than  in 
a  premixed  flame,  a  result  of  the  strong  correlation  between  thickening  and 
diffusion  as  expected.  Both  species  show  an  increase  in  total  production  rate 
due  to  increased  diffusion,  but  still  not  comparable  to  the  rate  of  production 
on  a  much  rehned  grid.  Generally  speaking,  the  thickened  flame  model  is 
favorable  to  combustion.  Although  applying  the  model  tends  to  postpone 
the  ignition  time,  it  allows  for  more  diffusion,  thus  more  combustion  without 
the  needs  to  resort  to  a  more  rehned  grid  for  hame  resolution. 

As  the  discussion  implies,  the  standard  used  in  this  context  to  determine 
the  accuracy  of  thickened  hame  and  reduced  mechanism  model  are  diherent 
for  premixed  and  dihusion  hames  in  certain  respects.  The  two  types  of  hames 
have  very  diherent  structures  and  combustion  dynamics.  Reactions  inside  a 
premixed  hame  regime  turn  out  to  be  more  volatile,  mainly  attributed  to  the 
interactions  (through  mass  dihusion  and  heat  transfer)  between  the  reaction 


zone  and  the  zone  with  burnt  gas  (discussed  in  Section  5.5.2).  The  hame 
speed  is  generally  considered  to  be  a  more  suitable  parameter  for  character¬ 
izing  the  dynamics  of  a  premixed  hame.  In  the  case  of  a  dihusion  hame, 
on  the  other  hand,  the  case  is  reversed.  While  the  hame  speed  is  not  well- 
dehned,  the  heat  release  is  a  quantihable  parameter  because  of  the  restricted 
region  the  combustion  could  take  place.  Detailed  analysis  shows  that  the 
reaction  rates  and  heat  release  are  both  negligible  outside  of  the  dihusion 
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Table  5.6:  The  le 
for  the  combinatic: 
tests. 


if  ignition  delay 
'  parameters  th 


^f^s  the  percentage  of  the  total  time) 
passed  the  premixed  flame  resolntion 


Full  chemical  mechanism 

Ax 

0.01  cm 

0.02  cm 

0.025  cm 

0.03  cm 

0.04  cm 

0.05  cm 

TF 

2 

3 

5 

6 

6 

8 

Ignition  delay 

7% 

15% 

29% 

36% 

36% 

50% 

Reduced  chemical  mechanism 

Ax 

0.01  cm 

0.02  cm 

0.025  cm 

0.03  cm 

0.04  cm 

0.05  cm 

TF 

1 

1 

2 

2 

2 

4 

Ignition  delay 

0 

0 

15% 

16% 

16% 

43% 

flame  reaction  zone  (the  length  of  which  is  dehned  as  the  flame  width). 

In  snmmary,  while  it  is  shown  in  the  case  of  a  premixed  flame  that  a 
too  small  thickening  factor  may  fail  to  resolve  the  flame,  onr  stndies  farther 
indicates  that  for  large  TF,  the  strnctnre  of  the  flame  and  the  combustion 
can  also  be  effected  by  delayed  burning.  Combining  the  cases  investigated 
above,  a  grid  size  of  Ax  =  0.05  cm  should  still  be  considered  as  too  coarse 
for  resolving  a  flame  that  combines  both  premixed  and  non-premixed  flame 
structures.  Using  a.  TF  =  8  would  underpredict  the  level  of  combustion, 
a  case  to  be  further  reviewed  in  Section  [6l  With  one  more  level  of  mesh 
rehnement,  a  grid  size  of  Ax  =  0.025  cm  is  suggested  to  give  a  more  or  less 
satisfactory  prediction  of  the  flame  structure  and  burning  as  a  whole,  with  a 


thickening  factor  of  TF  =  2  for  a  reduced  chemistry  mechanism.  Table  5.6 


further  shows  the  ignition  delay  for  the  whole  range  of  grid  sizes  investigated 
in  the  premixed  flame  section. 
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5. 5. 3. 3  The  Effect  of  the  Reduced  Chemical  Mechanism  in  One-di¬ 


mensional  Framework  A  detailed  analysis  of  the  rate  of  each  reaction  in 
the  fnll  chemical  mechanism  shows  that  each  of  these  reactions  reaches  their 
maximnm  valne  aronnd  the  ignition  time,  where  the  temperatnre  (or  heat 
release)  prohle  achieves  largest  gradient.  Among  them,  the  hydrogen-oxygen 
chain  reactions  ^(^3,  ^2  and  in  the  fnll  mechanism  Table  5.1|are  the  ones 
with  the  highest  rates.  After  ignition  time  is  reached,  the  dominant  reactions 
in  terms  of  reaction  rates  are  ^(^0,  #1,  ^2,  and  ^3,  while  all  other  reactions 
have  rates  mnch  smaller,  in  the  range  of  one  to  three  magnitnde  of  order. 
In  terms  of  heat  release,  the  combination  step  reaction  ^6  contribntes  the 
largest  amonnt  of  thermal  energy  in  forming  the  major  prodnct  species  H2O. 

The  large  difference  in  the  range  of  reaction  rates  jnstihes  a  redaction 
of  the  detailed  mechanism.  As  discussed  in  Section  5.2.2[  globally  reduced 
mechanisms  have  at  most  six  steps  in  total  in  the  case  of  hydrogen-oxygen 
flames.  Various  types  and  different  levels  of  reduced  chemical  mechanisms 
can  be  found  in  literature,  for  example  in  10  T5, 29  ,  ^46  .  These  developed 
mechanisms  are  generally  fit  to  be  in  agreement  with  either  experimental 
data  or  simulation  results  with  detailed  mechanism  for  very  different  flame 
structures  and  combustion  processes  ranging  from  fuel-lean  to  fuel-rich  con- 
hgurations,  low  to  high  temperatures,  deflagration,  detonation  or  flame  ex¬ 
tinction  processes,  normal  (latm)  to  high  (over  10  atm)  pressure,  and  so  on. 
In  this  study,  a  few  of  the  four-step  reduced  chemical  mechanisms  are  tested, 
and  the  one  found  to  be  most  consistent  with  the  configuration  studied  here 


is  the  one  given  in  the  reduced  mechanism  Table  5.3 


The  investigated  four-step  reduced  mechanism  consists  of  three  hydrogen- 
oxygen  chain  reaction  steps  (#1,  ^^2,  ^^3  in  the  detailed  mechanism  table) 
and  one  recombination  step  (#5  in  the  table).  After  a  linear  combination 
of  the  elementary  intermediate  steps,  the  rates  of  the  resulting  reactions 
are  comparable  and  have  about  the  same  order  of  magnitude,  among  which 
reaction  in  reduced  mechanism  produces  the  largest  amount  of  heat. 
The  previous  discussion  of  a  premixed  flame  in  Section  |5.3.1  shows  that 
the  reduced  mechanism  generally  provides  satisfactory  agreement  in  terms 
of  predicting  flame  speed  and  major  characteristics  of  the  flame.  However, 


as  Figure  [5i4|  suggests,  ignition  time  inconsistency  can  still  be  seen  with  the 
reduced  mechanism,  when  the  same  level  of  thickening  is  compared.  The 
general  approach  of  establishing  a  reduced  mechanism  is  associated  with 
leaving  out  reactions  that  have  minor  contributions  to  the  whole  scheme, 
or  combining  fast  reactions  into  major  reaction  steps.  During  the  reduction 
process,  the  accuracy  in  heat  release  and  related  effects  are  therefore  reduced 
as  well. 

The  major  advantage  of  various  reduced  mechanisms  is  their  efficiency 
in  computational  cost.  Indeed,  no  matter  what  model  is  chosen,  flamelet  or 
finite  rate  chemistry,  the  combustion  process  modeling  is  alway  an  expensive 
procedure.  With  the  flame  configuration  discussed  in  this  chapter,  a  more 
than  10%  reduction  in  running  one  flow  through  in  Scramjet  time  scale  is 
observed. 
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5. 5. 3. 4  Conclusion  This  section  discusses  the  hnite  rate  chemistry  ap¬ 
proach  being  applied  to  the  Scramjet  conhguration  that  is  reduced  to  the 
one-dimensional  case.  The  approach  is  evaluated  together  with  the  use  of 
the  thickened  flame  model  that  artihcially  expands  the  flame  front,  and  with 
a  reduced  chemical  mechanism  that  also  reduces  the  time  complexity  of  the 
whole  simulation.  It  is  shown  here  that  the  thickened  flame  model  effec¬ 
tively  increases  the  flame  width  while  retaining  the  major  flame  features 
with  enough  accuracy.  This  allows  for  the  flame  to  be  resolved  on  a  rel¬ 
atively  coarse  grid  size.  In  another  direction,  the  feasibility  of  a  specihc 
four-step  reduced  chemical  mechanism  is  investigated.  The  result  indicates 
that  the  reduced  chemical  mechanism  generally  agrees  quite  well  with  de¬ 
tailed  mechanism,  with  slight  delay  in  predicted  ignition  time.  However,  a 
combination  of  the  mechanism  with  thickened  flame  model  generally  dete¬ 
riorates  the  simulation  results.  The  level  of  allowed  thickening  in  order  to 
retain  the  same  amount  of  accuracy  will  be  further  restricted.  In  this  case, 
more  levels  of  mesh  rehnement  will  be  needed  in  order  to  resolve  the  flame. 
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Figure  6.1:  Sketch  of  a  scramjet  combustion  chamber 

Chapter  6 


6  Scramjet  Modeling 


6.1  Problem  Configurations 

The  Scramjet  studied  in  this  work  is  an  experimental  Mach  7  aircraft. 
Combustion  occurs  within  a  combustion  chamber  with  the  flow  still  super¬ 
sonic,  but  at  a  reduced  Mach  number  M  =  2.4.  A  model  scramjet  combustor 
has  been  developed,  and  experiments  conducted  in  a  Stanford  6”  Expansion 
Tube  Facility.  The  flow  parameters  are  chosen  to  be  representative  of  the 
flow  conditions  found  in  a  combustion  chamber.  The  details  of  the  experi¬ 


ment  setup  can  be  found  in  45  .  Results  from  experiments  and  numerical 


simulations  conducted  at  Stanford  University  Predictive  Science  Academic 
Alliance  Program  (PSAAP)  are  compared  with  simulation  results  presented 
in  this  chapter. 

A  sketch  of  a  simplihed  Scramjet  combustion  chamber  is  shown  in  Figure 
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6.1  The  dimensions  of  the  model  combustion  chamber  is  31cm  x  7.5  cm  x 


2.5  cm.  There  is  an  injection  port  with  a  diameter  of  0.2  cm  located  at  the 
center  line  of  the  top  wall,  with  a  distance  to  the  flat  plate  leading  edge  L  = 
7.06  cm.  The  H2  flows  out  of  the  nozzle  vertically  and  is  bent  downstream  by 
the  crossflow  O2  stream.  The  lower  wall  of  the  model  combustor  consists  of  a 
ramp  section  with  an  angle  of  10°.  Six  high-bandwidth  pressure  transducers 
are  mounted  in  plugs  inserted  into  the  centerline  of  the  lower  wall  of  the 


model  combustor  to  allow  pressure  measurements  45 


The  inflowing  oxygen  has  an  initial  temperature  Tq^  =  1548  K,  pressure 
P02  =  0.4  bar  and  crossflow  density  U02  =  180  cm/ms.  Together  with  a 
kinematic  viscosity  oi  v  =  5.36  x  10“^  cm^/ms,  the  Reynolds  number  for  the 
crossflow  can  be  computed  as  Reo2  =  LU jv  =  2.37  x  10^,  and  a  Kolmogorov 
scale  rjK,02  =  LRe~^^‘^  =  7.3  microns. 

The  hydrogen  is  injected  into  the  combustion  chamber  with  a  pressure 
Ph2  =  12.5  bar,  temperature  Th2  =  300  K  and  jet  exit  velocity  Uh2  = 
113.25  cm/ms.  Mixing  occurs  at  the  edges  of  the  plume,  leading  to  com¬ 
bustion  of  2H2  -|-  O2  — H2O  through  the  detailed  chemistry  mechanism  with 
eight  reacting  species  and  a  total  of  42  reactions.  With  the  kinematic  vis¬ 
cosity  u  =  1.6  X  10“^cm^/ms,  the  Reynolds  number  =  1-4  x  10®  and 

corresponding  Kolmogorov  scale  'r]K,H2  =  0-27 micron  are  calculated. 

The  jet-to-crossflow  momentum  ratio  dehned  as 


J  = 


fuel 


PoxidizerU^, 


(6.1) 


oxidizer 
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characterizes  the  general  structure  of  a  jet  in  cross-flow  conhguration,  and  is 
a  parameter  that  directly  affect  the  extent  of  penetration.  Experiments  such 
as  has  been  conducted  to  investigate  the  detailed  flow  structures  with 
varying  values  of  J.  To  be  compatible  with  the  setup  given  there,  J  =  5  is 
used  in  our  simulations. 

We  use  a  grid  spacing  of  0.1  cm  x  0.05  cm  x  0.05  cm  in  this  simulation.  It  is 
designed  to  be  more  rehned  on  y  and  2;  directions  to  reflect  the  fact  that  the 
mean  velocity  in  x  direction  dominates  the  direction  of  the  flow,  therefore 
the  characteristic  length  scale  is  much  larger  in  the  x  than  in  the  y  or  z 
directions.  To  increase  the  level  of  inflowing  turbulence  and  the  combustion 
that  follows,  a  randomized  turbulence  is  imposed  on  the  injected  crossflow, 
with  a  2.3%  intensity. 

To  model  the  flow  turbulent  motions  and  complicated  shock  structures, 
we  use  a  hfth-order  WENO  scheme  based  on  the  hnite  volume  discretiza¬ 


tion  with  an  explicit  third-order  TVD  Runge-Kutta  scheme  (Section  4.2), 


together  with  the  dynamic  LES  approach  (Section  3.3[).  The  combustion  is 


characterized  by  the  hnite  rate  chemistry  model  (Section  5.4.2). 


6.2  Results  and  Discussions 

Within  the  chamber,  the  turbulent  hame  is  neither  perfectly  premixed 
nor  perfectly  non-premixed,  displaying  the  complexities  of  pure  mixing,  par¬ 
tially  and  perfectly  mixed  combustion.  Based  on  the  analysis  in  Section 


5.5.3,  a  thickening  parameter  of  TF  =  8  is  chosen  for  the  current  grid  size. 
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It  has  already  been  shown  there  that  with  the  finite  rate  chemistry  model, 
the  current  grid  resolution  applied  to  three-dimensional  Scramjet  combus¬ 
tion  (Ao:  =  0.05  cm)  would  still  under-resolve  some  major  flame  structures, 
including  the  amount  of  heat  release,  the  start  of  ignition,  and  major  species 
production,  like  OH.  Here,  we  compare  the  simulation  results  to  the  ones 
conducted  with  the  same  parameter  configurations  at  the  Center  for  Turbu¬ 
lence  Research  at  Stanford  University.  The  Stanford  simulation  is  based  on 
a  fiamelet  description  of  the  flame  structure,  and  uses  an  energy-conserving 
unstructured  grid.  Therefore,  under  comparison  here  are  two  approaches  to 
modeling  turbulent  combustion,  namely,  the  finite  rate  chemistry  approach 


and  the  fiamelet  model.  Figure  |6.2|  and  Figure  |6.3|  compares  the  OH  mass 
fraction  produced  using  the  two  models.  Both  the  x  —  y  and  x  —  z  cross- 
sectional  profiles  show  that  combustion  starts  and  develops  in  the  chamber, 
though  not  as  well  developed  as  the  result  of  the  Stanford  simulation. 


Q 


The  major  discrepancy  lies  in  the  limited  burning  area,  especially  the 
first  half  of  the  injected  hydrogen  stream.  This  phenomena  is  largely  due  to 


the  delay  in  ignitid^jThe  one- dimensional  analysis  predicts  that  a  thicken¬ 


ing  factor|-^j^  =  8  would  delay  the  burning  time  by  about  half  of  the  total 


flow  time.  This  is  consistent  with  both  cross-sectional  profiles  of  the  OH 


production.  Figure  6.4  and  Figure  6.5  further  illustrate  the  contour  of  tem¬ 
perature  on  the  same  two  cross-sectional  planes,  and  further  show  the  sign 
of  insufficient  burning  in  the  first  half  of  the  hydrogen  downstream  flow. 


Figure  6.6  shows  the  comparison  of  the  lower  wall  pressure  among  exper- 
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OH; 


0.01  0.03  0.05  0.07  0.09  0.11  0.13  0.15 


Figure  6.2:  The  OH  mass  fraction  profiles  in  an  x  —  7/  plane,  the  slice  is 
taken  at  the  middle  of  the  height  of  the  combustor.  Comparison  is  between 
the  finite  rate  model  (upper)  and  the  flamelet  model  (lower). 


OH:  0.01  0.03  0.05  0.07  0.09  0.11  0.13 


Figure  6.3:  The  OH  mass  fraction  profile  in  an  x  —  ^  plane,  the  slice  is  taken 
at  the  middle  of  the  depth  of  the  combustor.  Comparison  is  between  the 
hnite  rate  model  (upper)  and  the  flamelet  model  (lower). 


imental  data  (measurements  taken  at  ten  locations  on  the  wall  by  the  pres¬ 
sure  transducers),  the  finite  rate  chemistry  model  and  the  flamelet  model  by 
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Figure  6.4:  The  temperature  profile  in  an  a;  —  ?/  plane,  the  slice  is  taken  at 
the  middle  of  the  height  of  the  combustor.  Comparison  is  between  the  finite 
rate  model  (upper)  and  the  flamelet  model  (lower). 


Figure  6.5:  The  temperature  profile  in  an  a;  —  2;  plane,  the  slice  is  taken  at 
the  middle  of  the  depth  of  the  combustor.  Comparison  is  between  the  hnite 
rate  model  (upper)  and  the  flamelet  model  (lower). 


Stanford,  evaluated  on  two  grid  levels.  It  shows  that  except  for  the  initiih 
combustion  period  (13cm  <  x  <  18cm),  the  finite  rate  model  agrees  with 
the  flamelet  model  and  experimental  data  in  most  places,  and  is  more  inline 
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with  the  coarse  grid  level  of  flamelet  model  up  to  around  13  cm.  However, 

>0.  of  .es- 


as  stated  above,  the  delayed  ignition  time  leads  to  about  odto 
idence  time  in  the  combustion  chamber  for  the  flame  portion  of  the  flow. 
In  addition,  due  to  the  nonlinear  coupling  between  the  chemistry  and  the 
hydrodynamics,  the  temperature  and  the  velocities  prohles  are  not  predicted 


correctly.  Further  ai^j^sis  shows  that  the  velocities  are  over  predicted.  The 


stream  wise  velocity  in  particular,  is  about  two  times  higher  than  predicted 
by  Stanford’s  simulation  results.  This  leads  to 


irther  decrease  in  the  com¬ 


bustion  and  heat  release  as  well.  On  the  other  hand,  the  boundary  layer 
drag  on  the  walls  is  about  double  the  amount  observed  in  the  Stanford’s  cal¬ 
culation.  The  combinational  effect  therefore  leads  to  the  agreement  of  wall 
pressures  with  experimental  data,  but  due  to  canceling  errors. 


Table  5T  indicates  the  minimum  mesh  requirements  after  applying  a  com¬ 
bination  of  the  thickened  flame  model  and  a  possible  reduced  chemistry  level. 
The  delay  in  ignition  could  be  reduced  to  7%  and  0%  for  a  full  and  reduced 
chemical  mechanism,  respectively,  within  the  considered  range  of  mesh  sizes. 
Based  on  the  one  dimensional  analysis,  at  least  the  second  case  would  lead 
to  a  correct  calculation.  The  predicted  minimum  mesh  size  without  a  de¬ 
lay  in  ignition  is  therefore  a  Ax  =  0.02  cm  with  reduced  chemistry,  and  a 
Ax  =  0.01  cm  with  full  chemistry  and  thickening  factor  TF  =  1  (no  thicken¬ 
ing). 
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Lower  Wall  Pressure  Profile 


Figure  6.6:  The  lower  wall  pressure,  compared  among  experimental  data, 
finite  rate  chemistry  model  and  the  fiamelet  model  by  Stanford  evaluated  on 
two  grid  levels. 


Chapter  7 


7  Conclusions 

In  this  work,  we  study  the  turbulent  mixing  and  combustion  of  the  su¬ 
personic  flow  using  Large  Eddy  Simulation.  Our  main  focus  is  to  determine 
mesh  requirements  needed  to  employ  finite  rate  chemistry  model.  The  com- 
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bustion  process  modeled  by  the  finite  rate  chemistry  approach  is  based  on 
the  notion  that  chemical  processes  have  length  scales  orders  of  magnitude 
larger  than  the  smallest  turbulence  scales,  and  therefore  could  be  resolved 

con- 


on  a  coarser  grid.  The  one-dimensional  studies  presented  in  Section  5.5.3 


firms  that  chemistry  could  be  resolved  with  a  grid  on  the  order  of  10“^  cm. 
Moreover,  to  compensate  for  the  deficiencies  in  resolving  combustion,  and  to 
reduce  computational  cost  (a  challenge  to  be  faced  when  using  more  refined 
grids),  it  could  be  useful  to  apply  the  thickened  flame  model  and  a  reduced 
chemical  mechanism.  Chapter  6  presents  a  three-dimensional  Scramjet  sim¬ 
ulation  with  the  grid  resolution  Ax  =  0.05  cm.  It  can  be  concluded  from 
the  simulation  results  that  while  correctly  predicting  the  major  features  of 
hydrodynamics  and  pressure  levels  at  the  nozzle  exit  (thrust).  Ax  =  0.05  cm 
is  still  not  refined  sufficiently  to  obtain  a  fully  developed  flame.  To  reach  that 
objective,  our  one-dimensional  flame  analysis  predicts  that  one  more  level  of 
mesh  refinement  to  Ax  =  0.025  cm  with  a  combination  of  thickened  flame 
model  could  succeed. 


Q 
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